From c30a0c45dedba0549018b5b0839f0e7b3966d783 Mon Sep 17 00:00:00 2001
From: Andrea-Havron-NOAA
 <85530309+Andrea-Havron-NOAA@users.noreply.github.com>
Date: Mon, 26 Aug 2024 20:40:59 +0000
Subject: [PATCH] modify test files after main change

---
 .../testthat/helper-integration-tests-setup.R |   3 +-
 .../test-integration-fims-estimation.R        | 938 ++++++------------
 2 files changed, 288 insertions(+), 653 deletions(-)

diff --git a/tests/testthat/helper-integration-tests-setup.R b/tests/testthat/helper-integration-tests-setup.R
index 531a51eb..50c87fa2 100644
--- a/tests/testthat/helper-integration-tests-setup.R
+++ b/tests/testthat/helper-integration-tests-setup.R
@@ -78,7 +78,7 @@ setup_and_run_FIMS <- function(iter_id,
   # recruit deviations should enter the model in normal space.
   # The log is taken in the likelihood calculations
   # alternative setting: recruitment$log_devs <- rep(0, length(om_input$logR.resid))
-  recruitment$log_devs <- methods::new(ParameterVector, om_input$logR.resid, om_input$nyr)
+  recruitment$log_devs <- methods::new(ParameterVector, om_input$logR.resid[-1], om_input$nyr-1)
 
   recruitment_distribution <- new(TMBDnormDistribution)
   # set up logR_sd using the normal log_sd parameter
@@ -254,6 +254,7 @@ setup_and_run_FIMS <- function(iter_id,
   return(list(
     parameters = parameters,
     obj = obj,
+    opt = opt,
     report = report,
     sdr_report = sdr_report,
     sdr_fixed = sdr_fixed,
diff --git a/tests/testthat/test-integration-fims-estimation.R b/tests/testthat/test-integration-fims-estimation.R
index cb85d9e6..c9bcf7b8 100644
--- a/tests/testthat/test-integration-fims-estimation.R
+++ b/tests/testthat/test-integration-fims-estimation.R
@@ -1,287 +1,221 @@
 load(test_path("fixtures", "integration_test_data.RData"))
 
-# test_that("deterministic test of fims", {
-#   # run function defined above to set up the test environment
-#   deterministic_env <- setup_fims(
-#     om_input = om_input,
-#     om_output = om_output,
-#     em_input = em_input
-#   )
-#   # Set-up
-#   deterministic_env$fims$CreateTMBModel()
-#   # CreateTMBModel calls a function in information that loops
-#   # over all the populations and fleets and sets all the pointers
-#   parameters <- list(p = deterministic_env$fims$get_fixed())
-#   # get_fixed function is an Rcpp function that loops over all Rcpp
-#   # modules and returned a vector of parameters being estimated
-
-#   # Set up TMB's computational graph
-#   obj <- MakeADFun(data = list(), parameters, DLL = "FIMS")
-
-#   # Calculate standard errors
-#   sdr <- TMB::sdreport(obj)
-#   sdr_fixed <- summary(sdr, "fixed")
-
-#   # Call report using deterministic parameter values
-#   # obj$report() requires parameter list to avoid errors
-#   report <- obj$report(obj$par)
-
-#   # Compare log(R0) to true value
-#   fims_logR0 <- sdr_fixed[1, "Estimate"]
-#   expect_gt(fims_logR0, 0.0)
-#   expect_equal(fims_logR0, log(om_input$R0))
-
-#   # Compare numbers at age to true value
-#   for (i in 1:length(c(t(om_output$N.age)))) {
-#     expect_equal(report$naa[[1]][i], c(t(om_output$N.age))[i])
-#   }
-
-#   # Compare biomass to true value
-#   for (i in 1:length(om_output$biomass.mt)) {
-#     expect_equal(report$biomass[[1]][i], om_output$biomass.mt[i])
-#   }
-
-#   # Compare spawning biomass to true value
-#   for (i in 1:length(om_output$SSB)) {
-#     expect_equal(report$ssb[[1]][i], om_output$SSB[i])
-#   }
-
-#   # Compare recruitment to true value
-#   fims_naa <- matrix(report$naa[[1]][1:(om_input$nyr * om_input$nages)],
-#     nrow = om_input$nyr, byrow = TRUE
-#   )
-
-#   # loop over years to compare recruitment by year
-#   for (i in 1:om_input$nyr) {
-#     expect_equal(fims_naa[i, 1], om_output$N.age[i, 1])
-#   }
-
-#   # confirm that recruitment matches the numbers in the first age
-#   # by comparing to fims_naa (what's reported from FIMS)
-#   expect_equal(
-#     fims_naa[1:om_input$nyr, 1],
-#     report$recruitment[[1]][1:om_input$nyr]
-#   )
-
-#   # confirm that recruitment matches the numbers in the first age
-#   # by comparing to the true values from the OM
-#   for (i in 1:om_input$nyr) {
-#     expect_equal(report$recruitment[[1]][i], om_output$N.age[i, 1])
-#   }
-
-#   # recruitment log_devs (fixed at initial "true" values)
-#   # the initial value of om_input$logR.resid is dropped from the model
-#   expect_equal(report$log_recruit_dev[[1]], om_input$logR.resid[-1])
-
-#   # F (fixed at initial "true" values)
-#   expect_equal(report$F_mort[[1]], om_output$f)
-
-#   # Expected catch
-#   fims_index <- report$exp_index
-#   for (i in 1:length(om_output$L.mt$fleet1)) {
-#     expect_equal(fims_index[[1]][i], om_output$L.mt$fleet1[i])
-#   }
-
-#   # Expect small relative error for deterministic test
-#   fims_object_are <- rep(0, length(em_input$L.obs$fleet1))
-#   for (i in 1:length(em_input$L.obs$fleet1)) {
-#     fims_object_are[i] <- abs(fims_index[[1]][i] - em_input$L.obs$fleet1[i]) / em_input$L.obs$fleet1[i]
-#   }
-
-#   # Expect 95% of relative error to be within 2*cv
-#   expect_lte(sum(fims_object_are > om_input$cv.L$fleet1 * 2.0), length(em_input$L.obs$fleet1) * 0.05)
-
-#   # Compare expected catch number at age to true values
-#   for (i in 1:length(c(t(om_output$L.age$fleet1)))) {
-#     expect_equal(report$cnaa[[1]][i], c(t(om_output$L.age$fleet1))[i])
-#   }
-
-#   # Expected catch number at age in proportion
-#   # QUESTION: Isn't this redundant with the non-proportion test above?
-#   fims_cnaa <- matrix(report$cnaa[[1]][1:(om_input$nyr * om_input$nages)],
-#     nrow = om_input$nyr, byrow = TRUE
-#   )
-#   fims_cnaa_proportion <- fims_cnaa / rowSums(fims_cnaa)
-#   om_cnaa_proportion <- om_output$L.age$fleet1 / rowSums(om_output$L.age$fleet1)
-
-#   for (i in 1:length(c(t(om_cnaa_proportion)))) {
-#     expect_equal(c(t(fims_cnaa_proportion))[i], c(t(om_cnaa_proportion))[i])
-#   }
-
-#   # Expected survey index.
-#   # Using [[2]] because the survey is the 2nd fleet.
-#   cwaa <- matrix(report$cwaa[[2]][1:(om_input$nyr * om_input$nages)], nrow = om_input$nyr, byrow = TRUE)
-#   expect_equal(fims_index[[2]], apply(cwaa, 1, sum) * om_output$survey_q$survey1)
-
-#   for (i in 1:length(om_output$survey_index_biomass$survey1)) {
-#     expect_equal(fims_index[[2]][i], om_output$survey_index_biomass$survey1[i])
-#   }
-
-#   fims_object_are <- rep(0, length(em_input$surveyB.obs$survey1))
-#   for (i in 1:length(em_input$survey.obs$survey1)) {
-#     fims_object_are[i] <- abs(fims_index[[2]][i] - em_input$surveyB.obs$survey1[i]) / em_input$surveyB.obs$survey1[i]
-#   }
-#   # Expect 95% of relative error to be within 2*cv
-#   expect_lte(sum(fims_object_are > om_input$cv.survey$survey1 * 2.0), length(em_input$surveyB.obs$survey1) * 0.05)
-
-#   # Expected catch number at age in proportion
-#   fims_cnaa <- matrix(report$cnaa[[2]][1:(om_input$nyr * om_input$nages)],
-#     nrow = om_input$nyr, byrow = TRUE
-#   )
-
-#   for (i in 1:length(c(t(om_output$survey_age_comp$survey1)))) {
-#     expect_equal(report$cnaa[[2]][i], c(t(om_output$survey_age_comp$survey1))[i])
-#   }
-
-#   fims_cnaa_proportion <- fims_cnaa / rowSums(fims_cnaa)
-#   om_cnaa_proportion <- om_output$survey_age_comp$survey1 / rowSums(om_output$survey_age_comp$survey1)
-
-#   for (i in 1:length(c(t(om_cnaa_proportion)))) {
-#     expect_equal(c(t(fims_cnaa_proportion))[i], c(t(om_cnaa_proportion))[i])
-#   }
-#   # clear memory
-#   deterministic_env$fims$clear()
-# })
-
-
-# test_that("nll test of fims", {
-#   nll_env <- setup_fims(
-#     om_input = om_input,
-#     om_output = om_output,
-#     em_input = em_input
-#   )
-#   # Set-up TMB
-#   nll_env$fims$CreateTMBModel()
-#   parameters <- list(p = nll_env$fims$get_fixed())
-#   par_list <- 1:length(parameters[[1]])
-#   par_list[2:length(par_list)] <- NA
-#   map <- list(p = factor(par_list))
-
-#   obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS", map = map)
-
-#   sdr <- TMB::sdreport(obj)
-#   sdr_fixed <- summary(sdr, "fixed")
-
-#   # log(R0)
-#   fims_logR0 <- sdr_fixed[1, "Estimate"]
-#   # expect_lte(abs(fims_logR0 - log(om_input$R0)) / log(om_input$R0), 0.0001)
-#   expect_equal(fims_logR0, log(om_input$R0))
-
-#   # Call report using deterministic parameter values
-#   # obj$report() requires parameter list to avoid errors
-#   report <- obj$report(obj$par)
-#   obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS", map = map)
-#   jnll <- obj$fn()
-
-
-#   # recruitment likelihood
-#   # log_devs is of length nyr-1
-#   rec_nll <- 0
-#   for(i in 1:(om_input$nyr - 1)){
-#     rec_nll <- rec_nll - dnorm(
-#       nll_env$recruitment$log_devs[i]$value, 0,
-#       om_input$logR_sd, TRUE
-#     )
-#   }
-
-#   # catch and survey index expected likelihoods
-#   index_nll_fleet <- -sum(dnorm(
-#     log(nll_env$catch),
-#     log(om_output$L.mt$fleet1),
-#     sqrt(log(em_input$cv.L$fleet1^2 + 1)), TRUE
-#   ))
-#   index_nll_survey <- -sum(dnorm(
-#     log(nll_env$survey_index),
-#     log(om_output$survey_index_biomass$survey1),
-#     sqrt(log(em_input$cv.survey$survey1^2 + 1)), TRUE
-#   ))
-#   index_nll <- index_nll_fleet + index_nll_survey
-#   # age comp likelihoods
-#   fishing_acomp_observed <- em_input$L.age.obs$fleet1
-#   fishing_acomp_expected <- om_output$L.age$fleet1 / rowSums(om_output$L.age$fleet1)
-#   survey_acomp_observed <- em_input$survey.age.obs$survey1
-#   survey_acomp_expected <- om_output$survey_age_comp$survey1 / rowSums(om_output$survey_age_comp$survey1)
-#   age_comp_nll_fleet <- age_comp_nll_survey <- 0
-#   for (y in 1:om_input$nyr) {
-#     age_comp_nll_fleet <- age_comp_nll_fleet -
-#       dmultinom(
-#         fishing_acomp_observed[y, ] * em_input$n.L$fleet1, em_input$n.L$fleet1,
-#         fishing_acomp_expected[y, ], TRUE
-#       )
-
-#     age_comp_nll_survey <- age_comp_nll_survey -
-#       dmultinom(
-#         survey_acomp_observed[y, ] * em_input$n.survey$survey1, em_input$n.survey$survey1,
-#         survey_acomp_expected[y, ], TRUE
-#       )
-#   }
-#   age_comp_nll <- age_comp_nll_fleet + age_comp_nll_survey
-#   expected_jnll <- rec_nll + index_nll + age_comp_nll
-
-#   expect_equal(report$nll_components[1], rec_nll)
-#   expect_equal(report$nll_components[2], index_nll_fleet)
-#   expect_equal(report$nll_components[3], age_comp_nll_fleet)
-#   expect_equal(report$nll_components[4], index_nll_survey)
-#   expect_equal(report$nll_components[5], age_comp_nll_survey)
-#   expect_equal(jnll, expected_jnll)
-
-#   nll_env$fims$clear()
-# })
-
-
-#   # recruitment likelihood
-#   # log_devs is of length nyr-1
-#   rec_nll <- 0
-#   for(i in 1:(om_input$nyr - 1)){
-#     rec_nll <- rec_nll - dnorm(
-#       nll_env$recruitment$log_devs[i]$value, 0,
-#       om_input$logR_sd, TRUE
-#     )
-#   }
-
-#   # catch and survey index expected likelihoods
-#   index_nll_fleet <- -sum(dnorm(
-#     log(nll_env$catch),
-#     log(om_output$L.mt$fleet1),
-#     sqrt(log(em_input$cv.L$fleet1^2 + 1)), TRUE
-#   ))
-#   index_nll_survey <- -sum(dnorm(
-#     log(nll_env$survey_index),
-#     log(om_output$survey_index_biomass$survey1),
-#     sqrt(log(em_input$cv.survey$survey1^2 + 1)), TRUE
-#   ))
-#   index_nll <- index_nll_fleet + index_nll_survey
-#   # age comp likelihoods
-#   fishing_acomp_observed <- em_input$L.age.obs$fleet1
-#   fishing_acomp_expected <- om_output$L.age$fleet1 / rowSums(om_output$L.age$fleet1)
-#   survey_acomp_observed <- em_input$survey.age.obs$survey1
-#   survey_acomp_expected <- om_output$survey_age_comp$survey1 / rowSums(om_output$survey_age_comp$survey1)
-#   age_comp_nll_fleet <- age_comp_nll_survey <- 0
-#   for (y in 1:om_input$nyr) {
-#     age_comp_nll_fleet <- age_comp_nll_fleet -
-#       dmultinom(
-#         fishing_acomp_observed[y, ] * em_input$n.L$fleet1, em_input$n.L$fleet1,
-#         fishing_acomp_expected[y, ], TRUE
-#       )
-
-#     age_comp_nll_survey <- age_comp_nll_survey -
-#       dmultinom(
-#         survey_acomp_observed[y, ] * em_input$n.survey$survey1, em_input$n.survey$survey1,
-#         survey_acomp_expected[y, ], TRUE
-#       )
-#   }
-#   age_comp_nll <- age_comp_nll_fleet + age_comp_nll_survey
-#   expected_jnll <- rec_nll + index_nll + age_comp_nll
-
-#   expect_equal(report$nll_components[1], rec_nll)
-#   expect_equal(report$nll_components[2], index_nll_fleet)
-#   expect_equal(report$nll_components[3], age_comp_nll_fleet)
-#   expect_equal(report$nll_components[4], index_nll_survey)
-#   expect_equal(report$nll_components[5], age_comp_nll_survey)
-#   expect_equal(jnll, expected_jnll)
-
-#   nll_env$fims$clear()
-# })
+test_that("deterministic test of fims", {
+  iter_id <- 1
+  result <- setup_and_run_FIMS(
+    iter_id = iter_id,
+    om_input_list = om_input_list,
+    om_output_list = om_output_list,
+    em_input_list = em_input_list,
+    estimation_mode = FALSE
+  )
+
+  # Set up TMB's computational graph
+  obj <- result$obj
+
+  # Calculate standard errors
+  sdr <- result$sdr
+  sdr_fixed <- result$sdr_fixed
+
+  # Call report using deterministic parameter values
+  # obj$report() requires parameter list to avoid errors
+  report <- result$report
+
+  # Compare log(R0) to true value
+  fims_logR0 <- sdr_fixed[1, "Estimate"]
+  expect_gt(fims_logR0, 0.0)
+  expect_equal(fims_logR0, log(om_input_list[[iter_id]]$R0))
+
+  # Compare numbers at age to true value
+  for (i in 1:length(c(t(om_output_list[[iter_id]]$N.age)))) {
+    expect_equal(report$naa[[1]][i], c(t(om_output_list[[iter_id]]$N.age))[i])
+  }
+
+  # Compare biomass to true value
+  for (i in 1:length(om_output_list[[iter_id]]$biomass.mt)) {
+    expect_equal(report$biomass[[1]][i], om_output_list[[iter_id]]$biomass.mt[i])
+  }
+
+  # Compare spawning biomass to true value
+  for (i in 1:length(om_output_list[[iter_id]]$SSB)) {
+    expect_equal(report$ssb[[1]][i], om_output_list[[iter_id]]$SSB[i])
+  }
+
+  # Compare recruitment to true value
+  fims_naa <- matrix(report$naa[[1]][1:(om_input_list[[iter_id]]$nyr * om_input_list[[iter_id]]$nages)],
+    nrow = om_input_list[[iter_id]]$nyr, byrow = TRUE
+  )
+
+  # loop over years to compare recruitment by year
+  for (i in 1:om_input_list[[iter_id]]$nyr) {
+    expect_equal(fims_naa[i, 1], om_output_list[[iter_id]]$N.age[i, 1])
+  }
+
+  # confirm that recruitment matches the numbers in the first age
+  # by comparing to fims_naa (what's reported from FIMS)
+  expect_equal(
+    fims_naa[1:om_input_list[[iter_id]]$nyr, 1],
+    report$recruitment[[1]][1:om_input_list[[iter_id]]$nyr]
+  )
+
+  # confirm that recruitment matches the numbers in the first age
+  # by comparing to the true values from the OM
+  for (i in 1:om_input_list[[iter_id]]$nyr) {
+    expect_equal(report$recruitment[[1]][i], om_output_list[[iter_id]]$N.age[i, 1])
+  }
+
+  # recruitment log_devs (fixed at initial "true" values)
+  # the initial value of om_input$logR.resid is dropped from the model
+  expect_equal(report$log_recruit_dev[[1]], om_input_list[[iter_id]]$logR.resid[-1])
+
+  # F (fixed at initial "true" values)
+  expect_equal(report$F_mort[[1]], om_output_list[[iter_id]]$f)
+
+  # Expected catch
+  fims_index <- report$exp_index
+  for (i in 1:length(om_output_list[[iter_id]]$L.mt$fleet1)) {
+    expect_equal(fims_index[[1]][i], om_output_list[[iter_id]]$L.mt$fleet1[i])
+  }
+
+  # Expect small relative error for deterministic test
+  fims_object_are <- rep(0, length(em_input_list[[iter_id]]$L.obs$fleet1))
+  for (i in 1:length(em_input_list[[iter_id]]$L.obs$fleet1)) {
+    fims_object_are[i] <- abs(fims_index[[1]][i] - em_input_list[[iter_id]]$L.obs$fleet1[i]) / em_input_list[[iter_id]]$L.obs$fleet1[i]
+  }
+
+  # Expect 95% of relative error to be within 2*cv
+  expect_lte(sum(fims_object_are > om_input_list[[iter_id]]$cv.L$fleet1 * 2.0), length(em_input_list[[iter_id]]$L.obs$fleet1) * 0.05)
+
+  # Compare expected catch number at age to true values
+  for (i in 1:length(c(t(om_output_list[[iter_id]]$L.age$fleet1)))) {
+    expect_equal(report$cnaa[[1]][i], c(t(om_output_list[[iter_id]]$L.age$fleet1))[i])
+  }
+
+  # Expected catch number at age in proportion
+  # QUESTION: Isn't this redundant with the non-proportion test above?
+  fims_cnaa <- matrix(report$cnaa[[1]][1:(om_input_list[[iter_id]]$nyr * om_input_list[[iter_id]]$nages)],
+    nrow = om_input_list[[iter_id]]$nyr, byrow = TRUE
+  )
+  fims_cnaa_proportion <- fims_cnaa / rowSums(fims_cnaa)
+  om_cnaa_proportion <- om_output_list[[iter_id]]$L.age$fleet1 / rowSums(om_output_list[[iter_id]]$L.age$fleet1)
+
+  for (i in 1:length(c(t(om_cnaa_proportion)))) {
+    expect_equal(c(t(fims_cnaa_proportion))[i], c(t(om_cnaa_proportion))[i])
+  }
+
+  # Expected survey index.
+  # Using [[2]] because the survey is the 2nd fleet.
+  cwaa <- matrix(report$cwaa[[2]][1:(om_input_list[[iter_id]]$nyr * om_input_list[[iter_id]]$nages)],
+    nrow = om_input_list[[iter_id]]$nyr, byrow = TRUE
+  )
+  expect_equal(fims_index[[2]], apply(cwaa, 1, sum) * om_output_list[[iter_id]]$survey_q$survey1)
+
+  for (i in 1:length(om_output_list[[iter_id]]$survey_index_biomass$survey1)) {
+    expect_equal(fims_index[[2]][i], om_output_list[[iter_id]]$survey_index_biomass$survey1[i])
+  }
+
+  fims_object_are <- rep(0, length(em_input_list[[iter_id]]$surveyB.obs$survey1))
+  for (i in 1:length(em_input_list[[iter_id]]$survey.obs$survey1)) {
+    fims_object_are[i] <- abs(fims_index[[2]][i] - em_input_list[[iter_id]]$surveyB.obs$survey1[i]) / em_input_list[[iter_id]]$surveyB.obs$survey1[i]
+  }
+  # Expect 95% of relative error to be within 2*cv
+  expect_lte(
+    sum(fims_object_are > om_input_list[[iter_id]]$cv.survey$survey1 * 2.0),
+    length(em_input_list[[iter_id]]$surveyB.obs$survey1) * 0.05
+  )
+
+  # Expected catch number at age in proportion
+  fims_cnaa <- matrix(report$cnaa[[2]][1:(om_input_list[[iter_id]]$nyr * om_input_list[[iter_id]]$nages)],
+    nrow = om_input_list[[iter_id]]$nyr, byrow = TRUE
+  )
+
+  for (i in 1:length(c(t(om_output_list[[iter_id]]$survey_age_comp$survey1)))) {
+    expect_equal(report$cnaa[[2]][i], c(t(om_output_list[[iter_id]]$survey_age_comp$survey1))[i])
+  }
+
+  fims_cnaa_proportion <- fims_cnaa / rowSums(fims_cnaa)
+  om_cnaa_proportion <- om_output_list[[iter_id]]$survey_age_comp$survey1 / rowSums(om_output_list[[iter_id]]$survey_age_comp$survey1)
+
+  for (i in 1:length(c(t(om_cnaa_proportion)))) {
+    expect_equal(c(t(fims_cnaa_proportion))[i], c(t(om_cnaa_proportion))[i])
+  }
+})
+
+
+test_that("nll test of fims", { 
+  iter_id <- 1
+
+  result <- setup_and_run_FIMS(
+    iter_id = iter_id,
+    om_input_list = om_input_list,
+    om_output_list = om_output_list,
+    em_input_list = em_input_list,
+    estimation_mode = FALSE
+  )
+
+  # Set up TMB's computational graph
+  obj <- result$obj
+  report <- result$report
+
+  # Calculate standard errors
+  sdr <- result$sdr
+  sdr_fixed <- result$sdr_fixed
+
+  # log(R0)
+  fims_logR0 <- sdr_fixed[1, "Estimate"]
+  # expect_lte(abs(fims_logR0 - log(om_input$R0)) / log(om_input$R0), 0.0001)
+  expect_equal(fims_logR0, log(om_input_list[[iter_id]]$R0))
+
+  # recruitment likelihood
+  # log_devs is of length nyr-1
+  rec_nll <- -sum(dnorm(
+    om_input_list[[iter_id]]$logR.resid[-1], rep(0, om_input_list[[iter_id]]$nyr-1),
+    om_input_list[[iter_id]]$logR_sd, TRUE
+  ))
+
+  # catch and survey index expected likelihoods
+  index_nll_fleet <- -sum(dlnorm(
+    em_input_list[[iter_id]]$L.obs$fleet1,
+    log(om_output_list[[iter_id]]$L.mt$fleet1),
+    sqrt(log(em_input_list[[iter_id]]$cv.L$fleet1^2 + 1)), TRUE
+  ))
+  index_nll_survey <- -sum(dlnorm(
+    em_input_list[[iter_id]]$surveyB.obs$survey1,
+    log(om_output_list[[iter_id]]$survey_index_biomass$survey1),
+    sqrt(log(em_input_list[[iter_id]]$cv.survey$survey1^2 + 1)), TRUE
+  ))
+  index_nll <- index_nll_fleet + index_nll_survey
+  # age comp likelihoods
+  fishing_acomp_observed <- em_input_list[[iter_id]]$L.age.obs$fleet1
+  fishing_acomp_expected <- om_output_list[[iter_id]]$L.age$fleet1 / rowSums(om_output_list[[iter_id]]$L.age$fleet1)
+  survey_acomp_observed <- em_input_list[[iter_id]]$survey.age.obs$survey1
+  survey_acomp_expected <- om_output_list[[iter_id]]$survey_age_comp$survey1 / rowSums(om_output_list[[iter_id]]$survey_age_comp$survey1)
+  age_comp_nll_fleet <- age_comp_nll_survey <- 0
+  for (y in 1:om_input_list[[iter_id]]$nyr) {
+    age_comp_nll_fleet <- age_comp_nll_fleet -
+      dmultinom(
+        fishing_acomp_observed[y, ] * em_input_list[[iter_id]]$n.L$fleet1, em_input_list[[iter_id]]$n.L$fleet1,
+        fishing_acomp_expected[y, ], TRUE
+      )
+
+    age_comp_nll_survey <- age_comp_nll_survey -
+      dmultinom(
+        survey_acomp_observed[y, ] * em_input_list[[iter_id]]$n.survey$survey1, em_input_list[[iter_id]]$n.survey$survey1,
+        survey_acomp_expected[y, ], TRUE
+      )
+  }
+  age_comp_nll <- age_comp_nll_fleet + age_comp_nll_survey
+  expected_jnll <- rec_nll + index_nll + age_comp_nll
+  jnll <- report$jnll
+
+  expect_equal(report$nll_components[1], rec_nll)
+  expect_equal(report$nll_components[2], index_nll_fleet)
+  expect_equal(report$nll_components[3], age_comp_nll_fleet)
+  expect_equal(report$nll_components[4], index_nll_survey)
+  expect_equal(report$nll_components[5], age_comp_nll_survey)
+  expect_equal(jnll, expected_jnll)
+
+})
 
 test_that("estimation test of fims", {
   # Initialize the iteration identifier and run FIMS with the 1st set of OM values
@@ -304,373 +238,73 @@ test_that("estimation test of fims", {
     em_input = em_input_list[[iter_id]]
   )
 })
-# test_that("run FIMS in a for loop with missing values", {
-#   for (i in 1:5) {
-#     fims <- Rcpp::Module("fims", PACKAGE = "FIMS")
-
-#     # Recruitment
-#     recruitment <- new(fims$BevertonHoltRecruitment)
-#     recruitment$log_rzero$value <- 13 # log(om_input$R0)
-#     # this change moves the starting value away from its true value
-#     recruitment$log_rzero$is_random_effect <- FALSE
-#     recruitment$log_rzero$estimated <- TRUE
-#     recruitment$logit_steep$value <- -log(1.0 - om_input$h) + log(om_input$h - 0.2)
-#     recruitment$logit_steep$is_random_effect <- FALSE
-#     recruitment$logit_steep$estimated <- FALSE
-#     recruitment$estimate_log_devs <- TRUE
-#     recruitment$log_devs <- methods::new(fims$ParameterVector, rep(0, length(om_input$logR.resid) - 1),length(om_input$logR.resid) - 1)
-
-#     recruitment_distribution <- new(fims$TMBDnormDistribution)
-# set dimension of observations
-#     recruitment_distribution$x <- new(ParameterVector, rep(0, om_input$nyr), om_input$nyr)
-#     recruitment_distribution$expected_values <- new(ParameterVector, rep(0, om_input$nyr), om_input$nyr)
-#     # set up logR_sd using the normal log_sd parameter
-#     # logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is
-#     # taken before the likelihood calculation
-#     recruitment_distribution$log_sd <- new(fims$ParameterVector, 1)
-#     recruitment_distribution$log_sd[1]$value <- log(om_input$logR_sd)
-#     recruitment_distribution$set_distribution_links("random_effects", recruitment$log_devs$get_id())
-
-#     # Data
-#     catch <- em_input$L.obs$fleet1
-#     fishing_fleet_index <- new(fims$Index, length(catch))
-#     fishing_fleet_index$index_data <- log(catch)
-#     testindex <- 2
-#     na_value <- -999
-#     if (i == 4) {
-#       fishing_fleet_index$index_data[testindex] <- na_value
-#     }
-#     fishing_fleet_age_comp <- new(fims$AgeComp, length(catch), om_input$nages)
-#     fishing_fleet_age_comp$age_comp_data <- c(t(em_input$L.age.obs$fleet1)) * em_input$n.L$fleet1
-#     if (i == 5) {
-#       fishing_fleet_age_comp$age_comp_data[testindex] <- na_value
-#     }
-
-#     survey_index <- em_input$surveyB.obs$survey1
-#     survey_fleet_index <- new(fims$Index, length(survey_index))
-#     survey_fleet_index$index_data <- log(survey_index)
-#     survey_fleet_age_comp <- new(fims$AgeComp, length(survey_index), om_input$nages)
-#     survey_fleet_age_comp$age_comp_data <- c(t(em_input$survey.age.obs$survey1)) * em_input$n.survey$survey1
-
-#     # Growth
-#     ewaa_growth <- new(fims$EWAAgrowth)
-#     ewaa_growth$ages <- om_input$ages
-#     ewaa_growth$weights <- om_input$W.mt
-
-#     # Maturity
-#     maturity <- new(fims$LogisticMaturity)
-#    maturity$inflection_point$value <- om_input$A50.mat
-#     maturity$inflection_point$is_random_effect <- FALSE
-#     maturity$inflection_point$estimated <- FALSE
-#     maturity$slope$value <- om_input$slope
-#     maturity$slope$is_random_effect <- FALSE
-#     maturity$slope$estimated <- FALSE
-
-#     # Fleet
-#     # Create the fishing fleet
-#     fishing_fleet_selectivity <- new(fims$LogisticSelectivity)
-#     fishing_fleet_selectivity$inflection_point$value <- om_input$sel_fleet$fleet1$A50.sel1
-#     fishing_fleet_selectivity$inflection_point$is_random_effect <- FALSE
-#     fishing_fleet_selectivity$inflection_point$estimated <- TRUE
-#     fishing_fleet_selectivity$slope$value <- om_input$sel_fleet$fleet1$slope.sel1
-#     fishing_fleet_selectivity$slope$is_random_effect <- FALSE
-#     fishing_fleet_selectivity$slope$estimated <- TRUE
-
-#     fishing_fleet <- new(fims$Fleet)
-#     fishing_fleet$nages <- om_input$nages
-#     fishing_fleet$nyears <- om_input$nyr
-#     fishing_fleet$log_Fmort <- new(fims$ParameterVector, log(om_output$f), om_input$nyr)
-#     fishing_fleet$log_Fmort$set_all_estimable(TRUE)
-#     #fishing_fleet$random_F <- FALSE
-#     fishing_fleet$log_q <- log(1.0)
-#     fishing_fleet$estimate_q <- FALSE
-#     fishing_fleet$random_q <- FALSE
-#     fishing_fleet$SetSelectivity(fishing_fleet_selectivity$get_id())
-
-#     # Set up fishery index data using the normal
-#     fishing_fleet_index_distribution <- methods::new(fims$TMBDnormDistribution)
-#     #normal observation error transformed on the log scale
-#     fishing_fleet_index_distribution$log_sd <- new(fims$ParameterVector, om_input$nyr)
-#     for(y in 1:om_input$nyr){
-#       fishing_fleet_index_distribution$log_sd[y]$value <- log(sqrt(log(em_input$cv.L$fleet1^2 + 1)))
-#     }
-#     fishing_fleet_index_distribution$log_sd$set_all_estimable(FALSE)
-#     # Set Data using the IDs from the modules defined above
-#     fishing_fleet_index_distribution$set_observed_data(fishing_fleet_index$get_id())
-#     fishing_fleet_index_distribution$set_distribution_links("data", fishing_fleet$log_expected_index$get_id())
-
-#     # Set up fishery age composition data using the multinomial
-#     fishing_fleet_agecomp_distribution <- methods::new(fims$TMBDmultinomDistribution)
-#     fishing_fleet_agecomp_distribution$set_observed_data(fishing_fleet_age_comp$get_id())
-#     fishing_fleet_agecomp_distribution$set_distribution_links("data", fishing_fleet$proportion_catch_numbers_at_age$get_id())
-
-
-#     # Create the survey fleet
-#     survey_fleet_selectivity <- new(fims$LogisticSelectivity)
-#     survey_fleet_selectivity$inflection_point$value <- om_input$sel_survey$survey1$A50.sel1
-#     survey_fleet_selectivity$inflection_point$is_random_effect <- FALSE
-#     survey_fleet_selectivity$inflection_point$estimated <- TRUE
-#     survey_fleet_selectivity$slope$value <- om_input$sel_survey$survey1$slope.sel1
-#     survey_fleet_selectivity$slope$is_random_effect <- FALSE
-#     survey_fleet_selectivity$slope$estimated <- TRUE
-
-#     survey_fleet <- new(fims$Fleet)
-#     survey_fleet$is_survey <- TRUE
-#     survey_fleet$nages <- om_input$nages
-#     survey_fleet$nyears <- om_input$nyr
-#     # survey_fleet$log_Fmort <- new(fims$ParameterVector, rep(log(0.0000000000000000000000000001), om_input$nyr), om_input$nyr) #-Inf?
-#     # survey_fleet$estimate_F <- FALSE
-#     # survey_fleet$random_F <- FALSE
-#     survey_fleet$log_q <- log(om_output$survey_q$survey1)
-#     survey_fleet$estimate_q <- TRUE
-#     survey_fleet$random_q <- FALSE
-#     survey_fleet$SetSelectivity(survey_fleet_selectivity$get_id())
-
-#     # Set up survey index data using the normal
-#     survey_fleet_index_distribution <- methods::new(fims$TMBDnormDistribution)
-#     #normal observation error transformed on the log scale
-#     # sd = sqrt(log(cv^2 + 1)), sd is log transformed
-#     survey_fleet_index_distribution$log_sd <- new(fims$ParameterVector, om_input$nyr)
-#     for(y in 1:om_input$nyr){
-#       survey_fleet_index_distribution$log_sd[y]$value <- log(sqrt(log(em_input$cv.survey$survey1^2 + 1)))
-#     }
-#     survey_fleet_index_distribution$log_sd$set_all_estimable(FALSE)
-#     # Set Data using the IDs from the modules defined above
-#     survey_fleet_index_distribution$set_observed_data(survey_fleet_index$get_id())
-#     survey_fleet_index_distribution$set_distribution_links("data", survey_fleet$log_expected_index$get_id())
-
-#     # Age composition data
-
-#     survey_fleet_agecomp_distribution <- methods::new(fims$TMBDmultinomDistribution)
-#     survey_fleet_agecomp_distribution$set_observed_data(survey_fleet_age_comp$get_id())
-#     survey_fleet_agecomp_distribution$set_distribution_links("data", survey_fleet$proportion_catch_numbers_at_age$get_id())
-
-
-#     # Population
-#     population <- new(fims$Population)
-#     # is it a problem these are not Parameters in the Population interface?
-#     # the Parameter class (from rcpp/rcpp_objects/rcpp_interface_base) cannot handle vectors,
-#     # do we need a ParameterVector class?
-#     population$log_M <- methods::new(fims$ParameterVector,
-#       rep(log(om_input$M.age[1]), om_input$nyr * om_input$nages),
-#       om_input$nyr * om_input$nages)
-#     population$log_M$set_all_estimable(FALSE)
-#     population$log_init_naa <- methods::new(fims$ParameterVector, log(om_output$N.age[1, ]), om_input$nages)
-#     population$log_init_naa$set_all_estimable(TRUE)
-#     population$nages <- om_input$nages
-#     population$ages <- om_input$ages
-#     population$nfleets <- sum(om_input$fleet_num, om_input$survey_num)
-#     population$nseasons <- 1
-#     population$nyears <- om_input$nyr
-#     population$SetMaturity(maturity$get_id())
-#     population$SetGrowth(ewaa_growth$get_id())
-#     population$SetRecruitment(recruitment$get_id())
-
-#     ## Set-up TMB
-#     fims$CreateTMBModel()
-#     parameters <- list(p = fims$get_fixed())
-#     obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS")
-
-#     opt <- with(obj, optim(par, fn, gr,
-#       method = "BFGS",
-#       control = list(maxit = 1000000, reltol = 1e-15)
-#     ))
-
-#     report <- obj$report(obj$par)
-#     g <- as.numeric(obj$gr(opt$par))
-#     h <- optimHess(opt$par, fn = obj$fn, gr = obj$gr)
-#     opt$par <- opt$par - solve(h, g)
-#     expect_false(is.null(report))
-
-#     max_gradient <- max(abs(obj$gr(opt$par)))
-#     expect_lte(max_gradient, 0.0001)
-#     fims$clear()
-#   }
-# })
-
-# test_that("agecomp in proportion works", {
-#   n.L_original <- om_input$n.L$fleet1
-#   n.survey_original <- om_input$n.survey$survey1
-#   om_input$n.L$fleet1 <- 1
-#   om_input$n.survey$survey1 <- 1
-#   on.exit(om_input$n.L$fleet1 <- n.L_original, add = TRUE)
-#   on.exit(om_input$n.survey$survey1 <- n.survey_original, add = TRUE)
-
-#   agecomp_in_proportion_env <- setup_fims(
-#     om_input = om_input,
-#     om_output = om_output,
-#     em_input = em_input
-#   )
-
-#   # Set-up TMB
-#   agecomp_in_proportion_env$fims$CreateTMBModel()
-#   # Create parameter list from Rcpp modules
-#   parameters <- list(p = agecomp_in_proportion_env$fims$get_fixed())
-#   obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS")
-
-#   opt <- with(obj, optim(par, fn, gr,
-#     method = "BFGS",
-#     control = list(maxit = 1000000, reltol = 1e-15)
-#   ))
-
-#   # Call report using MLE parameter values
-#   # obj$report() requires parameter list to avoid errors
-#   report <- obj$report(obj$env$last.par.best)
-#   sdr <- TMB::sdreport(obj)
-#   sdr_report <- summary(sdr, "report")
-#   # Numbers at age
-#   # Estimates and SE for NAA
-#   sdr_naa <- sdr_report[which(rownames(sdr_report) == "NAA"), ]
-#   naa_are <- rep(0, length(c(t(om_output$N.age))))
-#   for (i in 1:length(c(t(om_output$N.age)))) {
-#     naa_are[i] <- abs(sdr_naa[i, 1] - c(t(om_output$N.age))[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of NAA
-#   expect_lte(
-#     sum(naa_are > qnorm(.975) * sdr_naa[1:length(c(t(om_output$N.age))), 2]),
-#     0.05 * length(c(t(om_output$N.age)))
-#   )
-
-#   # Biomass
-#   sdr_biomass <- sdr_report[which(rownames(sdr_report) == "Biomass"), ]
-#   biomass_are <- rep(0, length(om_output$biomass.mt))
-#   for (i in 1:length(om_output$biomass.mt)) {
-#     biomass_are[i] <- abs(sdr_biomass[i, 1] - om_output$biomass.mt[i]) # / om_output$biomass.mt[i]
-#     # expect_lte(biomass_are[i], 0.15)
-#   }
-#   expect_lte(
-#     sum(biomass_are > qnorm(.975) * sdr_biomass[1:length(om_output$biomass.mt), 2]),
-#     0.05 * length(om_output$biomass.mt)
-#   )
-
-#   # Spawning biomass
-#   sdr_sb <- sdr_report[which(rownames(sdr_report) == "SSB"), ]
-#   sb_are <- rep(0, length(om_output$SSB))
-#   for (i in 1:length(om_output$SSB)) {
-#     sb_are[i] <- abs(sdr_sb[i, 1] - om_output$SSB[i]) # / om_output$SSB[i]
-#     # expect_lte(sb_are[i], 0.15)
-#   }
-#   expect_lte(
-#     sum(sb_are > qnorm(.975) * sdr_sb[1:length(om_output$SSB), 2]),
-#     0.05 * length(om_output$SSB)
-#   )
-
-#   # Recruitment
-#   fims_naa <- matrix(report$naa[[1]][1:(om_input$nyr * om_input$nages)],
-#     nrow = om_input$nyr, byrow = TRUE
-#   )
-#   sdr_naa1_vec <- sdr_report[which(rownames(sdr_report) == "NAA"), 2]
-#   sdr_naa1 <- sdr_naa1_vec[seq(1, om_input$nyr * om_input$nages, by = om_input$nages)]
-#   fims_naa1_are <- rep(0, om_input$nyr)
-#   for (i in 1:om_input$nyr) {
-#     fims_naa1_are[i] <- abs(fims_naa[i, 1] - om_output$N.age[i, 1]) # /
-#     # om_output$N.age[i, 1]
-#     # expect_lte(fims_naa1_are[i], 0.25)
-#   }
-#   expect_lte(
-#     sum(fims_naa1_are > qnorm(.975) * sdr_naa1[1:length(om_output$SSB)]),
-#     0.05 * length(om_output$SSB)
-#   )
-
-#   expect_equal(
-#     fims_naa[, 1],
-#     report$recruitment[[1]][1:om_input$nyr]
-#   )
-
-#   # recruitment log deviations
-#   # the initial value of om_input$logR.resid is dropped from the model
-#   sdr_rdev <- sdr_report[which(rownames(sdr_report) == "LogRecDev"), ]
-#   rdev_are <- rep(0, length(om_input$logR.resid) - 1)
-
-#   for (i in 1:(length(report$log_recruit_dev[[1]]) - 1)) {
-#     rdev_are[i] <- abs(report$log_recruit_dev[[1]][i] - om_input$logR.resid[i + 1]) # /
-#     #   exp(om_input$logR.resid[i])
-#     # expect_lte(rdev_are[i], 1) # 1
-#   }
-#   expect_lte(
-#     sum(rdev_are > qnorm(.975) * sdr_rdev[1:length(om_input$logR.resid) - 1, 2]),
-#     0.05 * length(om_input$logR.resid)
-#   )
-
-#   # F (needs to be updated when std.error is available)
-#   sdr_F <- sdr_report[which(rownames(sdr_report) == "FMort"), ]
-#   f_are <- rep(0, length(om_output$f))
-#   for (i in 1:length(om_output$f)) {
-#     f_are[i] <- abs(sdr_F[i, 1] - om_output$f[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of Fmort
-#   expect_lte(
-#     sum(f_are > qnorm(.975) * sdr_F[1:length(om_output$f), 2]),
-#     0.05 * length(om_output$f)
-#   )
-
-#   # Expected fishery catch and survey index
-#   fims_index <- sdr_report[which(rownames(sdr_report) == "ExpectedIndex"), ]
-#   fims_catch <- fims_index[1:om_input$nyr, ]
-#   fims_survey <- fims_index[(om_input$nyr + 1):(om_input$nyr * 2), ]
-
-#   # Expected fishery catch - om_output
-#   catch_are <- rep(0, length(om_output$L.mt$fleet1))
-#   for (i in 1:length(om_output$L.mt$fleet1)) {
-#     catch_are[i] <- abs(fims_catch[i, 1] - om_output$L.mt$fleet1[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of fishery catch
-#   expect_lte(
-#     sum(catch_are > qnorm(.975) * fims_catch[, 2]),
-#     0.05 * length(om_output$L.mt$fleet1)
-#   )
-
-#   # Expected fishery catch - em_input
-#   catch_are <- rep(0, length(em_input$L.obs$fleet1))
-#   for (i in 1:length(em_input$L.obs$fleet1)) {
-#     catch_are[i] <- abs(fims_catch[i, 1] - em_input$L.obs$fleet1[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of fishery catch
-#   expect_lte(
-#     sum(catch_are > qnorm(.975) * fims_catch[, 2]),
-#     0.05 * length(em_input$L.obs$fleet1)
-#   )
-
-
-#   # Expected fishery catch number at age
-#   sdr_cnaa <- sdr_report[which(rownames(sdr_report) == "CNAA"), ]
-#   cnaa_are <- rep(0, length(c(t(om_output$L.age$fleet1))))
-#   for (i in 1:length(c(t(om_output$L.age$fleet1)))) {
-#     cnaa_are[i] <- abs(sdr_cnaa[i, 1] - c(t(om_output$L.age$fleet1))[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of CNAA
-#   expect_lte(
-#     sum(cnaa_are > qnorm(.975) * sdr_cnaa[, 2]),
-#     0.05 * length(c(t(om_output$L.age$fleet1)))
-#   )
-#   print(cbind(sdr_cnaa, c(t(om_output$L.age$fleet1))))
-
-#   # Expected survey index - om_output
-#   index_are <- rep(0, length(om_output$survey_index_biomass$survey1))
-#   for (i in 1:length(om_output$survey_index_biomass$survey1)) {
-#     index_are[i] <- abs(fims_survey[i, 1] - om_output$survey_index_biomass$survey1[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of survey index
-#   expect_lte(
-#     sum(index_are > qnorm(.975) * fims_survey[, 2]),
-#     0.05 * length(om_output$survey_index_biomass$survey1)
-#   )
-
-#   # Expected survey index - em_input
-#   index_are <- rep(0, length(em_input$surveyB.obs$survey1))
-#   for (i in 1:length(em_input$surveyB.obs$survey1)) {
-#     index_are[i] <- abs(fims_survey[i, 1] - em_input$surveyB.obs$survey1[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of survey index
-#   # expect_lte(
-#   #   sum(index_are > qnorm(.975) * fims_survey[, 2]),
-#   #   0.05 * length(em_input$surveyB.obs$survey1)
-#   # )
-
-#   for (i in 1:length(em_input$surveyB.obs$survey1)) {
-#     expect_lte(abs(fims_survey[i, 1] - em_input$surveyB.obs$survey1[i]) /
-#       em_input$surveyB.obs$survey1[i], 0.25)
-#   }
-
-#   agecomp_in_proportion_env$fims$clear()
-# })
+
+test_that("run FIMS with missing values", {
+  # Initialize the iteration identifier
+  iter_id <- 1
+
+  # Define the NA (missing value) placeholder and the index where it will be inserted
+  na_value <- -999
+  na_index <- 2
+
+  # Introduce a missing value into the survey observations for the estimation model input
+  em_input_list[[iter_id]]$surveyB.obs$survey1[na_index] <- na_value
+
+  # Run the FIMS setup and execution function
+  result <- setup_and_run_FIMS(
+    iter_id = iter_id,
+    om_input_list = om_input_list,
+    om_output_list = om_output_list,
+    em_input_list = em_input_list,
+    estimation_mode = TRUE
+  )
+
+  # Validate that the result report is not null
+  expect_false(is.null(result$report))
+
+  # Obtain the gradient and Hessian matrix
+  g <- as.numeric(result$obj$gr(result$opt$par))
+  h <- optimHess(result$opt$par, fn = result$obj$fn, gr = result$obj$gr)
+  result$opt$par <- result$opt$par - solve(h, g)
+
+  # Obtain the maximum absolute gradient to check convergence
+  # Ensure that the maximum gradient is less than or equal to
+  # the specified tolerance (0.0001)
+  max_gradient <- max(abs(result$obj$gr(result$opt$par)))
+  expect_lte(max_gradient, 0.0001)
+})
+
+test_that("agecomp in proportion works", {
+  # Initialize the iteration identifier
+  iter_id <- 1
+
+  # Store the original values of the number of landings observations and
+  # survey observations
+  n.L_original <- om_input_list[[iter_id]]$n.L$fleet1
+  n.survey_original <- om_input_list[[iter_id]]$n.survey$survey1
+
+  # Set the number of landings observations and survey observations to 1
+  om_input_list[[iter_id]]$n.L$fleet1 <- 1
+  om_input_list[[iter_id]]$n.survey$survey1 <- 1
+  on.exit(om_input_list[[iter_id]]$n.L$fleet1 <- n.L_original, add = TRUE)
+  on.exit(om_input_list[[iter_id]]$n.survey$survey1 <- n.survey_original, add = TRUE)
+
+  # Run the FIMS setup and execution function
+  result <- setup_and_run_FIMS(
+    iter_id = iter_id,
+    om_input_list = om_input_list,
+    om_output_list = om_output_list,
+    em_input_list = em_input_list,
+    estimation_mode = TRUE
+  )
+
+  # Compare FIMS results with model comparison project OM values
+  validate_fims(
+    report = result$report,
+    sdr = TMB::sdreport(result$obj),
+    sdr_report = result$sdr_report,
+    om_input = om_input_list[[iter_id]],
+    om_output = om_output_list[[iter_id]],
+    em_input = em_input_list[[iter_id]]
+  )
+})