Adding the ability to customise prior parameters to all models. (#57)

* Updated stan files to have priors as a data argument for all top-level parameters.

* Changed affine model parameterisation. Updated vB model prior names to specify that max size adjust SD only. Updated models function to include default values.

* Added type specification to data structures for prior parameters because I am a dingus who forgot the first time almost everywhere.

* Removed misplaced code in Canham stan code.

* Updated vignettes to address the prior structure and talk about how to change the default values.

* Updated here be dragons short to comply with the new prior format. Moved the pre-print PDF to inst/include because it turns out that building the package locally deletes the inst/doc directory so there's a reason it was in the gitignore file.

* Put inst/doc/ back in the gitignore.

* Corrected index for SD prior ind_beta_1 so it's not defaulting to 0.

* Patching errors in testing.

* Correcting vignettes to load hmde before trying to use it.

* Added the used prior parameters to the generated quantities in order to allow the user to check values.

* Why oh why won't my generated quantites update?

* Trying different initialisation/assignment structure.

* Updated testing to correctly count the prior check parameters.

* Changed the assign data structure to account for different configurations of provided priors. Updated testing to reflect this.

* Added some assignment of plots.

.gitignore

@@ -31,6 +31,7 @@ vignettes/*.pdf
 # knitr and R markdown default cache directories
 *_cache/
 /cache/
+inst/doc/
 
 # Temporary files created by R markdown
 *.utf8.md

R/hmde_assign_data.R

@@ -15,7 +15,15 @@ hmde_assign_data <- function(model_template, data = NULL,...){
   if(!is.null(data)){ # Use provided tibble
     user_fields <- names(data)
 
+    user_code <- rlang::enquos(..., .check_assign = TRUE)
+    additional_user_fields <- names(user_code)
+    # Evaluate the RHS of expressions (the values)
+    additional_data <- purrr::map(user_code,
+                       ~rlang::eval_tidy(.x, env = rlang::caller_env())
+    )
+
   } else { # Grab user expressions from individual list items and extract data
+    additional_user_fields <- NULL
     user_code <- rlang::enquos(..., .check_assign = TRUE)
     user_fields <- names(user_code)
     # Evaluate the RHS of expressions (the values)
@@ -44,15 +52,21 @@ hmde_assign_data <- function(model_template, data = NULL,...){
   }
 
   for(i in model_fields){ # Iterate through required fields and fill them
-    if(i %in% user_fields){
+    if(i %in% user_fields){ #Check if the user has supplied it in a tibble
       model_template <- purrr::list_modify(model_template, !!!data[i])
-    } else if(is.null(model_template[[i]])){ # allows for default values
+
+    } else if(!is.null(additional_user_fields)){
+      if(i %in% additional_user_fields){ #Check if the user supplied it directly
+        model_template <- purrr::list_modify(model_template, !!!additional_data[i])
+      }
+    }
+
+    if(is.null(model_template[[i]])){ #Catches default tibble transformations
       model_template[[i]] <- switch(
         i,
         n_obs = length(data$y_obs),
         n_ind = length(unique(data$ind_id)),
-        y_bar = mean(data$y_obs),
-        model = model_template$model
+        y_bar = mean(data$y_obs)
       )
     }
 

R/hmde_models.R

@@ -34,6 +34,8 @@ hmde_const_single_ind <- function(){
        y_obs = NULL,
        obs_index = NULL,
        time = NULL,
+       prior_pars_ind_beta = c(0, 2),
+       prior_pars_global_error_sigma = c(0, 2),
        model = "constant_single_ind")
 }
 
@@ -48,6 +50,9 @@ hmde_const_multi_ind <- function(){
        obs_index = NULL,
        time = NULL,
        ind_id = NULL,
+       prior_pars_pop_beta_mu = c(0,2),
+       prior_pars_pop_beta_sigma = c(0,2),
+       prior_pars_global_error_sigma = c(0,2),
        model = "constant_multi_ind")
 }
 
@@ -60,6 +65,10 @@ hmde_canham_single_ind <- function(){
        y_obs = NULL,
        obs_index = NULL,
        time = NULL,
+       prior_pars_ind_max_growth = c(0,2),
+       prior_pars_ind_size_at_max_growth = c(0,2),
+       prior_pars_ind_k = c(0,2),
+       prior_pars_global_error_sigma = c(0,2),
        model = "canham_single_ind")
 }
 
@@ -74,6 +83,13 @@ hmde_canham_multi_ind <- function(){
        obs_index = NULL,
        time = NULL,
        ind_id = NULL,
+       prior_pars_pop_max_growth_mean = c(0,2),
+       prior_pars_pop_max_growth_sd = c(0,2),
+       prior_pars_pop_size_at_max_growth_mean = c(0,2),
+       prior_pars_pop_size_at_max_growth_sd = c(0,2),
+       prior_pars_pop_k_mean = c(0,2),
+       prior_pars_pop_k_sd = c(0,2),
+       prior_pars_global_error_sigma = c(0,2),
        model = "canham_multi_ind")
 }
 
@@ -88,6 +104,9 @@ hmde_vb_single_ind <- function(){
        obs_index = NULL,
        time = NULL,
        y_bar = NULL,
+       prior_pars_ind_max_size_sd_only = 2,
+       prior_pars_ind_growth_rate = c(0,2),
+       prior_pars_global_error_sigma = c(0,2),
        model = "vb_single_ind")
 }
 
@@ -103,6 +122,11 @@ hmde_vb_multi_ind <- function(){
        time = NULL,
        ind_id = NULL,
        y_bar = NULL,
+       prior_pars_pop_max_size_mean_sd_only = 2,
+       prior_pars_pop_max_size_sd = c(0,2),
+       prior_pars_pop_growth_rate_mean = c(0,2),
+       prior_pars_pop_growth_rate_sd = c(0,2),
+       prior_pars_global_error_sigma = c(0,2),
        model = "vb_multi_ind")
 }
 
@@ -118,7 +142,7 @@ hmde_affine_single_ind <- function(){
        time = NULL,
        int_method = NULL,
        y_bar = NULL,
-       prior_means = c(1,1),
-       prior_sds = c(2,2),
+       prior_pars_ind_const = c(1,2),
+       prior_pars_ind_beta_1 = c(0,2),
        model = "affine_single_ind")
 }

inst/stan/affine_single_ind.stan

@@ -66,8 +66,8 @@ data {
   real time[n_obs];
   real y_bar;
   int int_method;
-  real prior_means[2]; //vector of means for beta parameter priors
-  real prior_sds[2]; //Vector of SDs for beta parameter priors
+  real prior_pars_ind_const[2]; //vector of pars for beta_c parameter priors
+  real prior_pars_ind_beta_1[2]; //Vector of pars for beta_1 parameter priors
 }
 
 // The parameters accepted by the model.
@@ -122,16 +122,19 @@ model {
 
   //Priors
   //Individual level
-  ind_const ~normal(prior_means[1], prior_sds[1]);
-  ind_beta_1 ~lognormal(log(prior_means[2]), prior_sds[2]);
+  ind_const ~normal(prior_pars_ind_const[1], prior_pars_ind_const[2]);
+  ind_beta_1 ~lognormal(prior_pars_ind_beta_1[1], prior_pars_ind_beta_1[2]);
 }
 
 generated quantities{
   real y_hat[n_obs];
   array[3] real pars;
   real ind_beta_0;
   vector[1] y_temp;
-  int vers = 1;
+
+  //Return the used prior parameters
+  real check_prior_pars_ind_const[2] = prior_pars_ind_const;
+  real check_prior_pars_ind_beta[2] = prior_pars_ind_beta_1;
 
   ind_beta_0 = ind_const + ind_beta_1*y_bar;
 

inst/stan/canham_multi_ind.stan

@@ -16,6 +16,13 @@ data {
   int obs_index[n_obs];
   real time[n_obs];
   int ind_id[n_obs];
+  real prior_pars_pop_max_growth_mean[2];
+  real prior_pars_pop_max_growth_sd[2];
+  real prior_pars_pop_size_at_max_growth_mean[2];
+  real prior_pars_pop_size_at_max_growth_sd[2];
+  real prior_pars_pop_k_mean[2];
+  real prior_pars_pop_k_sd[2];
+  real prior_pars_global_error_sigma[2];
 }
 
 // The parameters accepted by the model.
@@ -74,21 +81,35 @@ model {
                    pop_k_sd);
 
   //Species level
-  pop_max_growth_mean ~normal(0, 1);
-  pop_max_growth_sd ~cauchy(0, 1);
-  pop_size_at_max_growth_mean ~normal(0, 1);
-  pop_size_at_max_growth_sd ~cauchy(0, 1);
-  pop_k_mean ~normal(0, 1);
-  pop_k_sd ~cauchy(0, 1);
+  pop_max_growth_mean ~normal(prior_pars_pop_max_growth_mean[1],
+                              prior_pars_pop_max_growth_mean[2]);
+  pop_max_growth_sd ~cauchy(prior_pars_pop_max_growth_sd[1],
+                            prior_pars_pop_max_growth_sd[2]);
+  pop_size_at_max_growth_mean ~normal(prior_pars_pop_size_at_max_growth_mean[1],
+                                      prior_pars_pop_size_at_max_growth_mean[2]);
+  pop_size_at_max_growth_sd ~cauchy(prior_pars_pop_size_at_max_growth_sd[1],
+                                    prior_pars_pop_size_at_max_growth_sd[2]);
+  pop_k_mean ~normal(prior_pars_pop_k_mean[1], prior_pars_pop_k_mean[2]);
+  pop_k_sd ~cauchy(prior_pars_pop_k_sd[1], prior_pars_pop_k_sd[2]);
 
   //Global level
-  global_error_sigma ~cauchy(0, 2);
+  global_error_sigma ~cauchy(prior_pars_global_error_sigma[1],
+                             prior_pars_global_error_sigma[2]);
 }
 
 generated quantities{
   real y_hat[n_obs];
   vector[1] y_temp;
 
+  //Return the used prior parameters
+  real check_prior_pars_pop_max_growth_mean[2] = prior_pars_pop_max_growth_mean;
+  real check_prior_pars_pop_max_growth_sd[2] = prior_pars_pop_max_growth_sd;
+  real check_prior_pars_pop_size_at_max_growth_mean[2] = prior_pars_pop_size_at_max_growth_mean;
+  real check_prior_pars_pop_size_at_max_growth_sd[2] = prior_pars_pop_size_at_max_growth_sd;
+  real check_prior_pars_pop_k_mean[2] = prior_pars_pop_k_mean;
+  real check_prior_pars_pop_k_sd[2] = prior_pars_pop_k_sd;
+  real check_prior_pars_global_error_sigma[2] = prior_pars_global_error_sigma;
+
   for(i in 1:n_obs){
 
     if(obs_index[i]==1){//Fits the first size

inst/stan/canham_single_ind.stan

@@ -14,6 +14,10 @@ data {
   real y_obs[n_obs];
   int obs_index[n_obs];
   real time[n_obs];
+  real prior_pars_ind_max_growth[2];
+  real prior_pars_ind_size_at_max_growth[2];
+  real prior_pars_ind_k[2];
+  real prior_pars_global_error_sigma[2];
 }
 
 // The parameters accepted by the model.
@@ -53,18 +57,27 @@ model {
 
   //Priors
   //Individual level
-  ind_max_growth ~lognormal(0, 1);
-  ind_size_at_max_growth ~lognormal(0, 1);
-  ind_k ~lognormal(0, 1);
+  ind_max_growth ~lognormal(prior_pars_ind_max_growth[1],
+                            prior_pars_ind_max_growth[2]);
+  ind_size_at_max_growth ~lognormal(prior_pars_ind_size_at_max_growth[1],
+                                    prior_pars_ind_size_at_max_growth[2]);
+  ind_k ~lognormal(prior_pars_ind_k[1], prior_pars_ind_k[2]);
 
   //Global level
-  global_error_sigma ~cauchy(0, 2);
+  global_error_sigma ~cauchy(prior_pars_global_error_sigma[1],
+                             prior_pars_global_error_sigma[2]);
 }
 
 generated quantities{
   real y_hat[n_obs];
   vector[1] y_temp;
 
+  //Return the used prior parameters
+  real check_prior_pars_ind_max_growth[2] = prior_pars_ind_max_growth;
+  real check_prior_pars_ind_size_at_max_growth[2] = prior_pars_ind_size_at_max_growth;
+  real check_prior_pars_ind_k[2] = prior_pars_ind_k;
+  real check_prior_pars_global_error_sigma[2] = prior_pars_global_error_sigma;
+
   for(i in 1:n_obs){
 
     if(obs_index[i]==1){//Fits the first size

inst/stan/constant_multi_ind.stan

@@ -18,6 +18,9 @@ data {
   int obs_index[n_obs];
   real time[n_obs];
   int ind_id[n_obs];
+  real prior_pars_pop_beta_mu[2];
+  real prior_pars_pop_beta_sigma[2];
+  real prior_pars_global_error_sigma[2];
 }
 
 // The parameters accepted by the model.
@@ -60,17 +63,29 @@ model {
                     pop_beta_sigma);
 
   //Population level
-  pop_beta_mu ~ normal(0, 1);
-  pop_beta_sigma ~cauchy(0, 1);
+  pop_beta_mu ~ normal(prior_pars_pop_beta_mu[1],
+                       prior_pars_pop_beta_mu[2]);
+  pop_beta_sigma ~cauchy(prior_pars_pop_beta_sigma[1],
+                         prior_pars_pop_beta_sigma[2]);
 
   //Global level
-  global_error_sigma ~cauchy(0, 2);
+  global_error_sigma ~cauchy(prior_pars_global_error_sigma[1],
+                             prior_pars_global_error_sigma[2]);
 }
 
 // The output
 generated quantities {
   real y_hat[n_obs];
 
+  //Return the used prior parameters
+  real check_prior_pars_pop_beta_mu[2];
+  real check_prior_pars_pop_beta_sigma[2];
+  real check_prior_pars_global_error_sigma[2];
+
+  check_prior_pars_pop_beta_mu = prior_pars_pop_beta_mu;
+  check_prior_pars_pop_beta_sigma = prior_pars_pop_beta_sigma;
+  check_prior_pars_global_error_sigma = prior_pars_global_error_sigma;
+
   for(i in 1:n_obs){
 
     //Fits the first size

inst/stan/constant_single_ind.stan

@@ -16,6 +16,8 @@ data {
   real y_obs[n_obs];
   int obs_index[n_obs];
   real time[n_obs];
+  real prior_pars_ind_beta[2];
+  real prior_pars_global_error_sigma[2];
 }
 
 // The parameters accepted by the model.
@@ -49,16 +51,22 @@ model {
 
   //Priors
   //Individual level
-  ind_beta ~ lognormal(0.1, 1);
+  ind_beta ~ lognormal(prior_pars_ind_beta[1],
+                       prior_pars_ind_beta[2]);
 
   //Global level
-  global_error_sigma ~cauchy(0, 2);
+  global_error_sigma ~cauchy(prior_pars_global_error_sigma[1],
+                             prior_pars_global_error_sigma[2]);
 }
 
 // The output
 generated quantities {
   real y_hat[n_obs];
 
+  //Return the used prior parameters
+  real check_prior_pars_ind_beta[2] = prior_pars_ind_beta;
+  real check_prior_pars_global_error_sigma[2] = prior_pars_global_error_sigma;
+
   for(i in 1:n_obs){
     //Fits the first size
     if(obs_index[i]==1){

inst/stan/vb_multi_ind.stan

@@ -16,6 +16,11 @@ data {
   real time[n_obs];
   int ind_id[n_obs];
   real y_bar;
+  real prior_pars_pop_max_size_mean_sd_only;
+  real prior_pars_pop_max_size_sd[2];
+  real prior_pars_pop_growth_rate_mean[2];
+  real prior_pars_pop_growth_rate_sd[2];
+  real prior_pars_global_error_sigma[2];
 }
 
 // The parameters accepted by the model.
@@ -67,19 +72,32 @@ model {
   ind_max_size ~lognormal(pop_max_size_mean, pop_max_size_sd);
 
   //Population level
-  pop_growth_rate_mean ~normal(0, 2);
-  pop_growth_rate_sd ~cauchy(0, 2);
-  pop_max_size_mean ~normal(max(log(y_obs)), 2);
-  pop_max_size_sd ~cauchy(0, 2);
+  pop_max_size_mean ~normal(log(max(y_obs)),
+                            prior_pars_pop_max_size_mean_sd_only);
+  pop_max_size_sd ~cauchy(prior_pars_pop_max_size_sd[1],
+                          prior_pars_pop_max_size_sd[2]);
+  pop_growth_rate_mean ~normal(prior_pars_pop_growth_rate_mean[1],
+                               prior_pars_pop_growth_rate_mean[2]);
+  pop_growth_rate_sd ~cauchy(prior_pars_pop_growth_rate_sd[1],
+                             prior_pars_pop_growth_rate_sd[2]);
 
   //Global level
-  global_error_sigma ~cauchy(0, 2);
+  global_error_sigma ~cauchy(prior_pars_global_error_sigma[1],
+                             prior_pars_global_error_sigma[2]);
 }
 
 generated quantities{
   real y_hat[n_obs];
   array[3] real pars;
 
+  //Return the used prior parameters
+  real check_prior_pars_pop_max_size_mean_sd_only = prior_pars_pop_max_size_mean_sd_only;
+  real check_prior_pars_pop_max_size_mean_mean_max_obs = log(max(y_obs));
+  real check_prior_pars_pop_max_size_sd[2] = prior_pars_pop_max_size_sd;
+  real check_prior_pars_pop_growth_rate_mean[2] = prior_pars_pop_growth_rate_mean;
+  real check_prior_pars_pop_growth_rate_sd[2] = prior_pars_pop_growth_rate_sd;
+  real check_prior_pars_global_error_sigma[2] = prior_pars_global_error_sigma;
+
   for(i in 1:n_obs){
     // Initialise the parameters for the observation
     pars[1] = ind_max_size[ind_id[i]] - y_bar;