Merge pull request #452 from NOAA-FIMS/385-refactor-consistent-use-of…

…-bias-correction-or-bias-adjustment-across-fims

Remove bias correction/adjustment

inst/include/distributions/functors/tmb_distributions.hpp

@@ -77,9 +77,7 @@ struct Dlnorm : public DistributionsBase<T> {
   T x;       /*!< observation */
   T meanlog; /*!< mean of the distribution of log(x) */
   T sdlog;   /*!< standard deviation of the distribution of log(x) */
-  bool do_bias_correction =
-      false; /*!< whether or not to bias correct the distribution */
-
+
   Dlnorm() : DistributionsBase<T>() {}
 
   /**
@@ -94,11 +92,9 @@ struct Dlnorm : public DistributionsBase<T> {
   virtual const T evaluate(const bool& do_log) {
     T logx = log(x);
     T nll;
-    if (do_bias_correction) {
-      nll = dnorm(logx, meanlog - pow(sdlog, 2) / 2, sdlog, true) - logx;
-    } else {
-      nll = dnorm(logx, meanlog, sdlog, true) - logx;
-    }
+
+    nll = dnorm(logx, meanlog, sdlog, true) - logx;
+
     if (do_log) {
       return nll;
     } else {

inst/include/interface/rcpp/rcpp_interface.hpp

@@ -188,10 +188,6 @@ RCPP_MODULE(fims) {
       .field("estimate_deviations",
              &BevertonHoltRecruitmentInterface::estimate_deviations)
       .method("get_id", &BevertonHoltRecruitmentInterface::get_id)
-      .field("recruitment_bias_adj",
-             &BevertonHoltRecruitmentInterface::recruit_bias_adjustment)
-      .field("use_bias_correction",
-             &BevertonHoltRecruitmentInterface::use_bias_correction)
       .field("log_sigma_recruit",
              &BevertonHoltRecruitmentInterface::log_sigma_recruit)
       .method("evaluate", &BevertonHoltRecruitmentInterface::evaluate)
@@ -293,9 +289,7 @@ RCPP_MODULE(fims) {
       .method("evaluate", &DlnormDistributionsInterface::evaluate)
       .field("x", &DlnormDistributionsInterface::x)
       .field("meanlog", &DlnormDistributionsInterface::meanlog)
-      .field("sdlog", &DlnormDistributionsInterface::sdlog)
-      .field("do_bias_correction",
-             &DlnormDistributionsInterface::do_bias_correction);
+      .field("sdlog", &DlnormDistributionsInterface::sdlog);
 
   Rcpp::class_<DmultinomDistributionsInterface>("TMBDmultinomDistribution")
       .constructor()

inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp

@@ -31,8 +31,6 @@ class RecruitmentInterfaceBase : public FIMSRcppInterfaceBase {
   // deviations*/
   /// static bool constrain_deviations; /**< whether or not the rec devs are
   /// constrained*/
-  // static std::vector<double> rec_bias_adj; /**< a vector of bias adjustment
-  // values*/
 
   RecruitmentInterfaceBase() {
     this->id = RecruitmentInterfaceBase::id_g++;
@@ -72,12 +70,9 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase {
   Parameter logit_steep;       /**< steepness or the productivity of the stock*/
   Parameter log_rzero;         /**< recruitment at unfished biomass */
   Parameter log_sigma_recruit; /**< the log of the stock recruit deviations */
-  Rcpp::NumericVector recruit_bias_adjustment; /**<vector bias adjustment*/
   Rcpp::NumericVector deviations;              /**< recruitment deviations*/
   bool estimate_deviations =
       false; /**< boolean describing whether to estimate */
-  bool use_bias_correction =
-      false; /**< boolean describing whether to do bias correction */
 
   BevertonHoltRecruitmentInterface() : RecruitmentInterfaceBase() {}
 
@@ -105,18 +100,11 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase {
 
     NLL.log_sigma_recruit = this->log_sigma_recruit.value;
     NLL.recruit_deviations.resize(deviations.size());  // Vector from TMB
-    NLL.recruit_bias_adjustment.resize(
-        recruit_bias_adjustment.size());  // Vector from TMB
     for (int i = 0; i < deviations.size(); i++) {
       NLL.recruit_deviations[i] = deviations[i];
-      NLL.recruit_bias_adjustment[i] = recruit_bias_adjustment[i];
     }
     FIMS_LOG << "Rec devs being passed to C++ are " << deviations << std::endl;
 
-    Rcout << "Rec bias adj being passed to C++ are " << recruit_bias_adjustment
-          << std::endl;
-
-    NLL.use_recruit_bias_adjustment = this->use_bias_correction;
     NLL.estimate_recruit_deviations = this->estimate_deviations;
     return NLL.evaluate_nll();
   }
@@ -170,7 +158,6 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase {
       }
     }
 
-    b0->use_recruit_bias_adjustment = this->use_bias_correction;
     // add to Information
     d0->recruitment_models[b0->id] = b0;
 
@@ -220,7 +207,6 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase {
       }
     }
 
-    b1->use_recruit_bias_adjustment = this->use_bias_correction;
     // add to Information
     d1->recruitment_models[b1->id] = b1;
 
@@ -270,7 +256,6 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase {
       }
     }
 
-    b2->use_recruit_bias_adjustment = this->use_bias_correction;
     // add to Information
     d2->recruitment_models[b2->id] = b2;
 
@@ -320,7 +305,6 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase {
       }
     }
 
-    b3->use_recruit_bias_adjustment = this->use_bias_correction;
     // add to Information
     d3->recruitment_models[b3->id] = b3;
 

inst/include/interface/rcpp/rcpp_objects/rcpp_tmb_distribution.hpp

@@ -161,8 +161,6 @@ class DlnormDistributionsInterface : public DistributionsInterfaceBase {
   Parameter x;       /*!< observation */
   Parameter meanlog; /*!< mean of the distribution of log(x) */
   Parameter sdlog;   /*!< standard deviation of the distribution of log(x) */
-  bool do_bias_correction; /*!< true if the lognormal should be bias corrected,
-                              default FALSE */
 
   DlnormDistributionsInterface() : DistributionsInterfaceBase() {}
 
@@ -185,7 +183,6 @@ class DlnormDistributionsInterface : public DistributionsInterfaceBase {
     dlnorm.x = this->x.value;
     dlnorm.meanlog = this->meanlog.value;
     dlnorm.sdlog = this->sdlog.value;
-    dlnorm.do_bias_correction = this->do_bias_correction;
     return dlnorm.evaluate(do_log);
   }
 

inst/include/population_dynamics/population/population.hpp

@@ -535,7 +535,7 @@ struct Population : public FIMSObject<Type> {
     /*
       Sets derived vectors to zero
       Performs parameters transformations
-      Sets recruitment deviations to mean 0, and then adds bias adjustment.
+      Sets recruitment deviations to mean 0.
      */
     Prepare();
     /*

inst/include/population_dynamics/recruitment/functors/recruitment_base.hpp

@@ -35,18 +35,7 @@ struct RecruitmentBase : public FIMSObject<Type> {
       recruit_deviations;            /*!< A vector of recruitment deviations */
   bool constrain_deviations = false; /*!< A flag to indicate if recruitment
                                  deviations are summing to zero or not */
-  typename ModelTraits<Type>::DataVector
-      recruit_bias_adjustment; /*!< A vector of bias adj values
-       (incorporating sigma_recruit)*/
-  // Initially fixing bias adjustment (b_y in collobarative
-  // workflow specification) to 1.0.
-  // In the future, this would be set by the user.
-  typename ModelTraits<Type>::DataVector
-      recruit_bias_adjustment_fraction; /*!< A vector of bias adjustment
-                                         fractions (on the 0 to 1 range)*/
-  bool use_recruit_bias_adjustment =
-      true; /*!< A flag to indicate if recruitment deviations are bias adjusted
-             */
+
   Type log_sigma_recruit; /*!< Log standard deviation of log recruitment
                        deviations */
   Type log_rzero;         /*!< Log of unexploited recruitment.*/
@@ -62,13 +51,10 @@ struct RecruitmentBase : public FIMSObject<Type> {
   virtual ~RecruitmentBase() {}
 
   /**
-   * @brief Prepares the recruitment bias adjustment vector.
+   * @brief Prepares the recruitment deviations vector.
    *
    */
   void Prepare() {
-    this->recruit_bias_adjustment_fraction.resize(
-        this->recruit_deviations.size());
-    this->recruit_bias_adjustment.resize(this->recruit_deviations.size());
     this->PrepareConstrainedDeviations();
   }
 
@@ -97,9 +83,6 @@ struct RecruitmentBase : public FIMSObject<Type> {
       for (size_t i = 0; i < this->recruit_deviations.size(); i++) {
         dnorm.x = fims::log(this->recruit_deviations[i]);
         dnorm.mean = 0.0;
-        if (this->use_recruit_bias_adjustment) {
-          dnorm.mean -= this->recruit_bias_adjustment[i];
-        }
         nll -= dnorm.evaluate(true);
       }
 #endif
@@ -131,30 +114,6 @@ struct RecruitmentBase : public FIMSObject<Type> {
     }
   }
 
-  /** @brief Prepare recruitment bias adjustment.
-   *  Based on Methot & Taylor (2011).
-   */
-  void PrepareBiasAdjustment() {
-    Type recruit_bias_adjustment_size = this->recruit_bias_adjustment.size();
-
-    if (!this->use_recruit_bias_adjustment) {
-      for (size_t i = 0; i < recruit_bias_adjustment_size; i++) {
-        this->recruit_bias_adjustment_fraction[i] = 1.0;
-        this->recruit_bias_adjustment[i] = 0.0;
-      }
-    } else {
-      for (size_t i = 0; i < recruit_bias_adjustment_size; i++) {
-        // Initially fixing bias adjustment (b_y in collobarative
-        // workflow specification) to 1.0.
-        // In the future, this would be set by the user.
-        this->recruit_bias_adjustment_fraction[i] = 1.0;
-        this->recruit_bias_adjustment[i] =
-            0.5 * fims::exp(this->log_sigma_recruit) *
-            fims::exp(this->log_sigma_recruit) *
-            this->recruit_bias_adjustment_fraction[i];
-      }
-    }
-  }
 };
 
 template <class Type>

tests/gtest/test_population_dynamics_recruitment_base.cpp

@@ -33,52 +33,4 @@ namespace
     }
   }
 
-  TEST(recruitment_bias_adjustment, PrepareBiasAdjustment_works)
-  {
-    fims::SRBevertonHolt<double> recruit;
-    recruit.log_sigma_recruit = fims::log(0.7);
-
-    // Test if use_recruit_bias_adjustment = false works
-    recruit.use_recruit_bias_adjustment = false;
-    recruit.recruit_bias_adjustment = {0.0, 0.0, 0.0};
-    recruit.recruit_bias_adjustment_fraction = {1.0, 1.0, 1.0};
-    recruit.PrepareBiasAdjustment();
-
-    std::vector<double> expected_bias_adjustment_false = {0.0, 0.0, 0.0};
-    for (int i = 0; i < recruit.recruit_bias_adjustment.size(); i++)
-    {
-      EXPECT_EQ(recruit.recruit_bias_adjustment[i],
-                expected_bias_adjustment_false[i]);
-    }
-
-    // Test if use_recruit_bias_adjustment = true works
-    recruit.use_recruit_bias_adjustment = true;
-    recruit.recruit_bias_adjustment = {0.0, 0.0, 0.0};
-    recruit.recruit_bias_adjustment_fraction = {1.0, 1.0, 1.0};
-    recruit.PrepareBiasAdjustment();
-
-    // R code to generate true values: 0.5 * 0.7^2 * c(1.0, 1.0, 1.0)
-    // 0.245 0.245 0.245
-    std::vector<double> expected_bias_adjustment_true = {0.245, 0.245, 0.245};
-    for (int i = 0; i < recruit.recruit_bias_adjustment.size(); i++)
-    {
-      EXPECT_NEAR(recruit.recruit_bias_adjustment[i],
-                  expected_bias_adjustment_true[i], 0.001);
-    }
-
-    // Test if recruit_bias_adjustment_fraction will be converted to 1.0
-    // when initial input value is 2.0
-    recruit.recruit_bias_adjustment_fraction = {2.0, 2.0, 2.0};
-    recruit.PrepareBiasAdjustment();
-
-    // R code to generate true values: 0.5 * 0.7^2 * c(1.0, 1.0, 1.0)
-    // 0.245 0.245 0.245
-    std::vector<double> expected_bias_adjustment = {0.245, 0.245, 0.245};
-    for (int i = 0; i < recruit.recruit_bias_adjustment.size(); i++)
-    {
-      EXPECT_NEAR(recruit.recruit_bias_adjustment[i],
-                  expected_bias_adjustment[i], 0.001);
-    }
-  }
-
 }

tests/integration/integration_class.hpp

@@ -571,7 +571,6 @@ class IntegrationTest {
                     std::cout << "'logR.resid' not found.\n";
                 }
             }
-            rec->use_recruit_bias_adjustment = false;
             pop.recruitment = rec;
 
             // set maturity

tests/testthat/test-recruitment-interface.R

@@ -17,8 +17,6 @@ test_that("Recruitment input settings work as expected", {
   recruitment$logit_steep$estimated <- TRUE
   recruitment$log_rzero$value <- log(r0)
   recruitment$log_sigma_recruit$value <- log(0.7)
-  recruitment$recruitment_bias_adj <- rep(1.0, 3)
-  recruitment$use_bias_correction <- FALSE
 
   expect_equal(recruitment$get_id(), 1)
   expect_equal(recruitment$logit_steep$value, 0.78845736)
@@ -44,15 +42,5 @@ test_that("Recruitment input settings work as expected", {
   recruitment$estimate_deviations <- TRUE
   expect_equal(recruitment$evaluate_nll(), expected = expected_nll)
 
-  recruitment$use_bias_correction <- TRUE
-  recruitment$recruitment_bias_adj <- rep(0.245, 3)
-  expected_nll <- -sum(log(stats::dnorm(log(devs), -0.245, 0.7)))
-  expect_equal(recruitment$evaluate_nll(), expected = expected_nll)
-
-  recruitment$use_bias_correction <- TRUE
-  recruitment$recruitment_bias_adj <- c(0.245, 0.2, 0.1)
-  expected_nll <- -sum(log(stats::dnorm(log(devs), c(-0.245, -0.2, -0.1), 0.7)))
-  expect_equal(recruitment$evaluate_nll(), expected = expected_nll)
-
   fims$clear()
 })

tests/testthat/test-tmb-distributions.R

@@ -34,7 +34,6 @@ test_that("dlnorm", {
   dlnorm_$x$value <- y
   dlnorm_$meanlog$value <- 0
   dlnorm_$sdlog$value <- 1
-  dlnorm_$do_bias_correction <- FALSE
   # evaluate the density and compare with R
   expect_equal(dlnorm_$evaluate(TRUE), stats::dlnorm(y, 0, 1, TRUE))
   expect_equal(dlnorm_$evaluate(FALSE), stats::dlnorm(y, 0, 1, FALSE))