Update recruitment, priors

DESCRIPTION

@@ -6,7 +6,8 @@ Authors@R: c(
     comment = c(ORCID = "0000-0002-3658-4517")),
     person("Joe", "Thorley", , "[email protected]", role = "aut",
            comment = c(ORCID = "0000-0002-7683-4592")),
-    person("John", "Boulanger", role = "aut"),
+    person("John", "Boulanger", role = "aut",
+           comment = c(ORCID = "0000-0001-8222-1445")),
     person("Ayla", "Pearson", , "[email protected]", role = "aut"),
     person(given = "Environment and Climate Change Canada",
              role = "cph")
@@ -22,10 +23,12 @@ Depends:
 Imports: 
     bboudata,
     chk,
+    dplyr,
     generics,
     ggplot2,
     glue,
     lifecycle,
+    magrittr,
     mcmcderive,
     mcmcr,
     newdata,

NAMESPACE

@@ -65,6 +65,7 @@ export(bb_plot_year_recruitment)
 export(bb_plot_year_survival)
 export(bb_plot_year_trend_recruitment)
 export(bb_plot_year_trend_survival)
+export(bb_predict_calf_cow_ratio)
 export(bb_predict_growth)
 export(bb_predict_lambda)
 export(bb_predict_population_change)
@@ -95,9 +96,14 @@ export(samples)
 export(tidy)
 import(chk)
 import(glue)
+import(magrittr)
 import(nimble)
 importFrom(bboudata,bbd_chk_data_recruitment)
 importFrom(bboudata,bbd_chk_data_survival)
+importFrom(dplyr,arrange)
+importFrom(dplyr,bind_rows)
+importFrom(dplyr,left_join)
+importFrom(dplyr,tibble)
 importFrom(generics,augment)
 importFrom(generics,glance)
 importFrom(generics,tidy)

R/data-recruitment.R

@@ -12,6 +12,18 @@ data_clean_recruitment <- function(data, quiet = FALSE) {
 data_prep_recruitment <- function(data, year_start) {
   data$CowsBulls <- data$Cows + data$Bulls
   data$Year <- caribou_year(data$Year, data$Month, year_start = year_start)
+  data <-
+    data %>%
+    dplyr::group_by(Year) %>%
+    dplyr::summarize(
+      Cows = sum(.data$Cows),
+      CowsBulls = sum(.data$CowsBulls),
+      UnknownAdults = sum(.data$UnknownAdults),
+      Yearlings = sum(.data$Yearlings),
+      Calves = sum(.data$Calves),
+      PopulationName = dplyr::first(.data$PopulationName)
+    ) %>%
+    dplyr::ungroup()
   data$Annual <- factor(data$Year)
 
   data

R/data-survival.R

@@ -5,16 +5,28 @@ data_clean_survival <- function(data, quiet) {
   data
 }
 
-data_prep_survival <- function(data, include_uncertain_morts, year_start) {
+data_prep_survival <- function(data, include_uncertain_morts,
+                               year_start) {
+  data$Mortalities <- data$MortalitiesCertain
+  if (include_uncertain_morts) {
+    data$Mortalities <- data$Mortalities + data$MortalitiesUncertain
+  }
   data$Year <- caribou_year(data$Year, data$Month, year_start = year_start)
   data$Annual <- factor(data$Year)
+
   # leaves month but sets factor levels to be caribou month for model
   nmonth <- length(unique(data$Month))
   data$Month <- factor(data$Month, levels = month_levels(year_start, n = nmonth))
-  data$Mortalities <- data$MortalitiesCertain
-  if (include_uncertain_morts) {
-    data$Mortalities <- data$Mortalities + data$MortalitiesUncertain
-  }
+
+  data
+}
+
+data_adjust_intercept <- function(data) {
+  year_b0 <- year_intercept(data)
+  levels <- unique(data$Year)
+  levels <- levels[levels != year_b0]
+  levels <- c(year_b0, levels)
+  data$Annual <- factor(data$Year, levels = levels)
   data
 }
 

R/derived.R

@@ -34,12 +34,12 @@ derived_expr_recruitment <- function(fit, year) {
     lik <- extract_lik(fit)
   }
   paste0("for(i in 1:length(Annual)) {
-  log(prediction[i]) <- ", lik, "\n}")
+  logit(prediction[i]) <- ", lik, "\n}")
 }
 
 derived_expr_recruitment_trend <- function() {
   "for(i in 1:length(Annual)) {
-  log(prediction[i]) <- b0 + bYear * Year[i]\n}"
+  logit(prediction[i]) <- b0 + bYear * Year[i]\n}"
 }
 
 derived_expr_survival_trend <- function() {

R/fit-recruitment.R

@@ -1,14 +1,14 @@
 #' Fit Recruitment Model
 #'
 #' Fit heirarchical Bayesian recruitment model using Nimble.
-#' 
+#'
 #' If the number of years is > `min_random_year`, a fixed-effects model is fit.
 #' Otherwise, a mixed-effects model is fit with random intercept for each year.
-#' If `year_trend` is TRUE and the number of years is > `min_random_year`, the model 
-#' will be fit with year as a continuous effect (i.e. trend) and no fixed effect of year. 
-#' If `year_trend` is TRUE and the number of years is <= `min_random_year`, the model 
-#' will be fit with year as a continuous effect and a random intercept for each year. 
-#' 
+#' If `year_trend` is TRUE and the number of years is > `min_random_year`, the model
+#' will be fit with year as a continuous effect (i.e. trend) and no fixed effect of year.
+#' If `year_trend` is TRUE and the number of years is <= `min_random_year`, the model
+#' will be fit with year as a continuous effect and a random intercept for each year.
+#'
 #' The start month of the Caribou year can be adjusted with `year_start`.
 #'
 #' @inheritParams params
@@ -94,18 +94,18 @@ bb_fit_recruitment <- function(
 #' Fit Recruitment Model with Maximum Likelihood
 #'
 #' Fit recruitment model with Maximum Likelihood using Nimble Laplace Approximation.
-#' 
+#'
 #' If the number of years is > `min_random_year`, a fixed-effects model is fit.
 #' Otherwise, a mixed-effects model is fit with random intercept for each year.
-#' If `year_trend` is TRUE and the number of years is > `min_random_year`, the model 
-#' will be fit with year as a continuous effect (i.e. trend) and no fixed effect of year. 
-#' If `year_trend` is TRUE and the number of years is <= `min_random_year`, the model 
+#' If `year_trend` is TRUE and the number of years is > `min_random_year`, the model
+#' will be fit with year as a continuous effect (i.e. trend) and no fixed effect of year.
+#' If `year_trend` is TRUE and the number of years is <= `min_random_year`, the model
 #' will be fit with year as a continuous effect and a random intercept for each year.
 #'
 #' Year effect can be excluded with `exclude_year`. This can be useful if the ML model is failing to converge.
-#' 
+#'
 #' The start month of the Caribou year can be adjusted with `year_start`.
-#' 
+#'
 #' @inheritParams params
 #' @return A list of the Nimble model object and Maximum Likelihood output with estimates and standard errors on the transformed scale.
 #' @export
@@ -139,10 +139,9 @@ bb_fit_recruitment_ml <- function(
   data <- model_data_recruitment(data, year_start = year_start, quiet = quiet)
   year_random <- data$datal$nAnnual >= min_random_year
   if (!year_random && year_trend) {
-    if(!quiet) message_trend_fixed()
+    if (!quiet) message_trend_fixed()
   }
 
-  priors <- priors_recruitment()
   model <- model_recruitment(
     data = data$datal,
     year_random = year_random,
@@ -164,7 +163,7 @@ bb_fit_recruitment_ml <- function(
 
   convergence_fail <- ml_converge_fail(fit) || ml_se_fail(fit)
   if (convergence_fail) {
-    if(!quiet) message_convergence_fail()
+    if (!quiet) message_convergence_fail()
   }
 
   fit <- fit$result

R/fit-survival.R

@@ -11,16 +11,16 @@ ml_se_fail <- function(x) {
 #' Fit Survival Model
 #'
 #' Fits hierarchical Bayesian survival model using Nimble.
-#' 
+#'
 #' If the number of years is > `min_random_year`, a fixed-effects model is fit.
 #' Otherwise, a mixed-effects model is fit with random intercept for each year.
-#' If `year_trend` is TRUE and the number of years is > `min_random_year`, the model 
-#' will be fit with year as a continuous effect (i.e. trend) and no fixed effect of year. 
-#' If `year_trend` is TRUE and the number of years is <= `min_random_year`, the model 
+#' If `year_trend` is TRUE and the number of years is > `min_random_year`, the model
+#' will be fit with year as a continuous effect (i.e. trend) and no fixed effect of year.
+#' If `year_trend` is TRUE and the number of years is <= `min_random_year`, the model
 #' will be fit with year as a continuous effect and a random intercept for each year.
 #'
 #' The model is always fit with random intercept for each month.
-#' 
+#'
 #' The start month of the Caribou year can be adjusted with `year_start`.
 #'
 #' @inheritParams params
@@ -107,20 +107,20 @@ bb_fit_survival <- function(data,
 #' Fit Survival Model with Maximum Likelihood
 #'
 #' Fits hierarchical survival model with Maximum Likelihood using Nimble Laplace approximation.
-#' 
+#'
 #' If the number of years is > `min_random_year`, a fixed-effects model is fit.
 #' Otherwise, a mixed-effects model is fit with random intercept for each year.
-#' If `year_trend` is TRUE and the number of years is > `min_random_year`, the model 
-#' will be fit with year as a continuous effect (i.e. trend) and no fixed effect of year. 
-#' If `year_trend` is TRUE and the number of years is <= `min_random_year`, the model 
+#' If `year_trend` is TRUE and the number of years is > `min_random_year`, the model
+#' will be fit with year as a continuous effect (i.e. trend) and no fixed effect of year.
+#' If `year_trend` is TRUE and the number of years is <= `min_random_year`, the model
 #' will be fit with year as a continuous effect and a random intercept for each year.
 #'
 #' The model is always fit with random intercept for each month.
-#' 
+#'
 #' Year effect can be excluded with `exclude_year`. This can be useful if the ML model is failing to converge.
-#' 
+#'
 #' The start month of the Caribou year can be adjusted with `year_start`.
-#' 
+#'
 #' @inheritParams params
 #' @return A list of the Nimble model object and Maximum Likelihood output with estimates and standard errors on the transformed scale.
 #' @export
@@ -150,14 +150,22 @@ bb_fit_survival_ml <- function(data,
   chk_null_or(inits, vld = vld_named)
   chk_flag(quiet)
 
-  data <-
-    model_data_survival(data,
-      include_uncertain_morts = include_uncertain_morts,
-      year_start = year_start, quiet = quiet
-    )
-  year_random <- data$datal$nAnnual >= min_random_year
+  # special treatment of intercept for ML fixed
+  data <- data_clean_survival(data, quiet = quiet)
+  data <- data_prep_survival(data,
+    include_uncertain_morts = include_uncertain_morts,
+    year_start = year_start
+  )
+  year_random <- length(unique(data$Year)) >= min_random_year
+  if (!year_random) {
+    data <- data_adjust_intercept(data)
+  }
+
+  datal <- data_list_survival(data)
+  data <- list(datal = datal, data = data)
+
   if (!year_random && year_trend && !exclude_year) {
-    if(!quiet) message_trend_fixed()
+    if (!quiet) message_trend_fixed()
   }
 
   model <-
@@ -178,7 +186,7 @@ bb_fit_survival_ml <- function(data,
 
   convergence_fail <- ml_converge_fail(fit) || ml_se_fail(fit)
   if (convergence_fail) {
-    if(!quiet) message_convergence_fail()
+    if (!quiet) message_convergence_fail()
   }
 
   fit <- fit$result

R/model-recruitment.R

@@ -8,7 +8,7 @@ model_data_recruitment <- function(data, year_start = year_start, quiet) {
 #' Build Nimble recruitment model.
 #'
 #' This is for use by developers.
-#' 
+#'
 #' @inheritParams params
 #' @param demographic_stochasticity A flag indicating whether to include demographic_stochasticity in the recruitment model.
 #' @export
@@ -69,16 +69,16 @@ model_recruitment <-
       if (year_trend) {
         if (year_random) {
           for (i in 1:nObs) {
-            log(eRecruitment[i]) <- b0 + bAnnual[Annual[i]] + bYear * Year[i]
+            logit(eRecruitment[i]) <- b0 + bAnnual[Annual[i]] + bYear * Year[i]
           }
         } else {
           for (i in 1:nObs) {
-            log(eRecruitment[i]) <- b0 + bYear * Year[i]
+            logit(eRecruitment[i]) <- b0 + bYear * Year[i]
           }
         }
       } else {
         for (i in 1:nObs) {
-          log(eRecruitment[i]) <- b0 + bAnnual[Annual[i]]
+          logit(eRecruitment[i]) <- b0 + bAnnual[Annual[i]]
         }
       }
 
@@ -95,9 +95,9 @@ model_recruitment <-
         }
       } else {
         for (i in 1:nObs) {
-          FemaleYearlings[i] <- yearling_female_proportion * Yearlings[i]
+          FemaleYearlings[i] <- round(yearling_female_proportion * Yearlings[i])
           OtherAdultsFemales[i] <-
-            adult_female_proportion * UnknownAdults[i]
+            round(adult_female_proportion * UnknownAdults[i])
         }
       }
 
@@ -107,7 +107,7 @@ model_recruitment <-
           # in original model because max cannot be used with ML
           ((FemaleYearlings[i] + Cows[i] + OtherAdultsFemales[i]) < 1) +
           FemaleYearlings[i] + Cows[i] + OtherAdultsFemales[i]
-        Calves[i] ~ dpois(eRecruitment[i] * AdultsFemales[i])
+        Calves[i] ~ dbin(eRecruitment[i], AdultsFemales[i])
       }
     })
 
@@ -122,10 +122,10 @@ model_recruitment <-
       buildDerivs = TRUE,
       name = "bboumodel_recruitment"
     )
-    
+
     # reset to user
     nimbleOptions(verbose = verbose)
     nimbleOptions(enableDerivs = enable_derivs)
-    
+
     model
   }

R/model-survival.R

@@ -1,4 +1,7 @@
-model_data_survival <- function(data, include_uncertain_morts, year_start, quiet) {
+model_data_survival <- function(data,
+                                include_uncertain_morts,
+                                year_start,
+                                quiet) {
   data <- data_clean_survival(data, quiet = quiet)
   data <- data_prep_survival(data,
     include_uncertain_morts = include_uncertain_morts,
@@ -11,7 +14,7 @@ model_data_survival <- function(data, include_uncertain_morts, year_start, quiet
 #' Build Nimble survival model.
 #'
 #' This is for use by developers.
-#' 
+#'
 #' @inheritParams params
 #' @param build_derivs A flag indicating whether to build derivatives Laplace approximation.
 #' @export
@@ -94,16 +97,16 @@ model_survival <- function(data,
   enable_derivs <- nimbleOptions("enableDerivs")
   nimbleOptions(verbose = FALSE)
   nimbleOptions(enableDerivs = TRUE)
-  
+
   model <- nimbleModel(code,
     constants = constants,
     buildDerivs = build_derivs,
     name = "bboumodel_survival"
   )
-  
+
   # reset to user
   nimbleOptions(verbose = verbose)
   nimbleOptions(enableDerivs = enable_derivs)
-  
+
   model
 }

R/namespace.R

@@ -1,6 +1,7 @@
-#' @import chk glue nimble
+#' @import chk glue nimble magrittr
 #' @importFrom scales percent
 #' @importFrom rescale rescale
+#' @importFrom dplyr arrange tibble bind_rows left_join
 #' @importFrom purrr quietly map map_dbl map_lgl
 #' @importFrom stats setNames update median coef predict nobs plogis runif qnorm logLik
 #' @importFrom ggplot2 autoplot .data ggplot aes geom_pointrange geom_hline scale_x_continuous scale_x_discrete scale_y_continuous expand_limits xlab geom_line geom_ribbon

R/nimble.R

@@ -32,7 +32,7 @@ run_nimble_ml <- function(model, inits, prior_inits, quiet) {
   cmodel <- compileNimble(model, resetFunctions = TRUE)
   params <- setupMargNodes(model)$paramNodes
   chk_subset(names(inits), params, x_name = "Names in `inits`")
-  
+
   # need to deal with special case when estimating adult_female_proportion
   # nimble does not guess default param nodes correctly
   cows <- map_lgl(params, function(x) grepl("Cows", x))