Power model add (#13)

Added the power law model code and testing suite.

* Added power model code and set tests to run on this branch.

* Typo in vector of branches to run tests.

* Removed this branch from the test workflow.

* Updated power law model which required some minor changes to variable references in helper functions for testing and data generation.

* Fixed what I broke when I changed the test functions.

* Changed power law model to have a population-level single value for y_bar, including the data generation process. Currently fails to properly estimate coefficient parameter for single individual.

* Fix error, Use map functions

* Corrected error in power single individual test and changed test helper functions to use map rather than for loop.

---------

Co-authored-by: Daniel Falster <[email protected]>

DESCRIPTION

@@ -35,6 +35,8 @@ LinkingTo:
 SystemRequirements: GNU make
 Suggests: 
     knitr,
+    purrr,
+    dplyr,
     rmarkdown,
     testthat (>= 3.0.0),
     withr,

R/rmot_models.R

@@ -17,7 +17,9 @@ rmot_model <- function(model=NULL){
                    constant_single_ind = rmot_const_single_ind(),
                    constant_multi_ind = rmot_const_multi_ind(),
                    canham_single_ind = rmot_canham_single_ind(),
-                   canham_multi_ind = rmot_canham_multi_ind())
+                   canham_multi_ind = rmot_canham_multi_ind(),
+                   power_single_ind = rmot_power_single_ind(),
+                   power_multi_ind = rmot_power_multi_ind())
 
   class(output) <- "rmot_object"
 
@@ -93,3 +95,34 @@ rmot_canham_multi_ind <- function(){
        model = "canham_multi_ind")
 }
 
+#' Data configuration template for power law growth single individual model
+#' @keywords internal
+#' @noRd
+
+rmot_power_single_ind <- function(){
+  list(step_size = NULL,
+       n_obs = NULL,
+       y_obs = NULL,
+       obs_index = NULL,
+       time = NULL,
+       y_bar = NULL,
+       y_0_obs = NULL,
+       model = "power_single_ind")
+}
+
+#' Data configuration template for power law growth single species model
+#' @keywords internal
+#' @noRd
+
+rmot_power_multi_ind <- function(){
+  list(step_size = NULL,
+       n_obs = NULL,
+       n_ind = NULL,
+       y_obs = NULL,
+       obs_index = NULL,
+       time = NULL,
+       ind_id = NULL,
+       y_bar = NULL,
+       y_0_obs = NULL,
+       model = "power_multi_ind")
+}

R/rmot_run.R

@@ -20,7 +20,9 @@ rmot_run <- function(model_template, ...) {
          constant_single_ind = rstan::sampling(stanmodels$constant_single_ind, data = model_template, ...),
          constant_multi_ind = rstan::sampling(stanmodels$constant_multi_ind, data = model_template, ...),
          canham_single_ind = rstan::sampling(stanmodels$canham_single_ind, data = model_template, ...),
-         canham_multi_ind = rstan::sampling(stanmodels$canham_multi_ind, data = model_template, ...))
+         canham_multi_ind = rstan::sampling(stanmodels$canham_multi_ind, data = model_template, ...),
+         power_single_ind = rstan::sampling(stanmodels$power_single_ind, data = model_template, ...),
+         power_multi_ind = rstan::sampling(stanmodels$power_multi_ind, data = model_template, ...))
 
   return(out)
 }

R/stanmodels.R

@@ -1,14 +1,16 @@
 # Generated by rstantools.  Do not edit by hand.
 
 # names of stan models
-stanmodels <- c("canham_multi_ind", "canham_single_ind", "constant_multi_ind", "constant_single_ind", "linear")
+stanmodels <- c("canham_multi_ind", "canham_single_ind", "constant_multi_ind", "constant_single_ind", "linear", "power_multi_ind", "power_single_ind")
 
 # load each stan module
 Rcpp::loadModule("stan_fit4canham_multi_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4canham_single_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4constant_multi_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4constant_single_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4linear_mod", what = TRUE)
+Rcpp::loadModule("stan_fit4power_multi_ind_mod", what = TRUE)
+Rcpp::loadModule("stan_fit4power_single_ind_mod", what = TRUE)
 
 # instantiate each stanmodel object
 stanmodels <- sapply(stanmodels, function(model_name) {

inst/stan/power_multi_ind.stan

@@ -0,0 +1,151 @@
+// Power law model with multiple individuals
+
+functions{
+  //DE function for use with Runge-Kutta method
+  //pars = [ind_coeff, ind_power]
+  real DE(real y, vector pars){
+    return pars[1] * pow((y/pars[3]), -pars[2]);
+  }
+
+  real rk4_step(real y, vector pars, real interval){
+    real k1;
+    real k2;
+    real k3;
+    real k4;
+    real y_hat;
+
+    k1 = DE(y, pars);
+    k2 = DE(y+interval*k1/2.0, pars);
+    k3 = DE(y+interval*k2/2.0, pars);
+    k4 = DE(y+interval*k3, pars);
+
+    y_hat = y + (1.0/6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4) * interval;
+
+    return y_hat;
+  }
+
+  real rk4(real y, vector pars, real interval, real step_size){
+    int steps;
+    real duration;
+    real y_hat;
+    real step_size_temp;
+
+    duration = 0;
+    y_hat = y;
+
+    while(duration < interval){
+      //Determine the relevant step size
+      step_size_temp = min([step_size, interval-duration]);
+
+      //Get next size estimate
+      y_hat = rk4_step(y_hat, pars, step_size_temp);
+
+      //Increment observed duration
+      duration = duration + step_size_temp;
+    }
+
+    return y_hat;
+  }
+}
+
+// Data structure
+data {
+  real step_size;
+  int n_obs;
+  int n_ind;
+  real y_obs[n_obs];
+  int obs_index[n_obs];
+  real time[n_obs];
+  int ind_id[n_obs];
+  real y_bar;
+  real y_0_obs[n_ind];
+}
+
+// The parameters accepted by the model.
+parameters {
+  //Individual level
+  real<lower=0> ind_y_0[n_ind];
+  real<lower=0> ind_coeff[n_ind];
+  real<lower=0> ind_power[n_ind];
+
+  //Species level
+  real species_coeff_mean;
+  real<lower=0> species_coeff_sd;
+  real species_power_mean;
+  real<lower=0> species_power_sd;
+
+  //Global level
+  real<lower=0> global_error_sigma;
+}
+
+// The model to be estimated.
+model {
+  real y_hat[n_obs];
+  vector[3] pars;
+
+  for(i in 1:n_obs){
+    // Initialise the parameters for the observation
+    pars[1] = ind_coeff[ind_id[i]];
+    pars[2] = ind_power[ind_id[i]];
+    pars[3] = y_bar;
+
+    if(obs_index[i]==1){//Fits the first size
+      y_hat[i] = ind_y_0[ind_id[i]];
+    }
+
+    if(i < n_obs){
+      if(ind_id[i+1]==ind_id[i]){
+        //Estimate next size
+        y_hat[i+1] = rk4(y_hat[i], pars, (time[i+1] - time[i]), step_size);
+      }
+    }
+  }
+
+  //Likelihood
+  y_obs ~ normal(y_hat, global_error_sigma);
+
+  //Priors
+  //Individual level
+  ind_y_0 ~ normal(y_0_obs, global_error_sigma);
+  ind_coeff ~lognormal(species_coeff_mean, species_coeff_sd);
+  ind_power ~lognormal(species_power_mean, species_power_sd);
+
+  //Species level
+  species_coeff_mean ~normal(0, 2);
+  species_coeff_sd ~cauchy(0, 2);
+  species_power_mean ~normal(0, 2);
+  species_power_sd ~cauchy(0, 2);
+
+  //Global level
+  global_error_sigma ~cauchy(0, 5);
+}
+
+generated quantities{
+  real y_hat[n_obs];
+  real Delta_hat[n_obs];
+  vector[3] pars;
+
+  for(i in 1:n_obs){
+    // Initialise the parameters for the observation
+    pars[1] = ind_coeff[ind_id[i]];
+    pars[2] = ind_power[ind_id[i]];
+    pars[3] = y_bar;
+
+    if(obs_index[i]==1){//Fits the first size
+      y_hat[i] = ind_y_0[ind_id[i]];
+    }
+
+    if(i < n_obs){
+      if(ind_id[i+1]==ind_id[i]){
+        //Estimate next size
+        y_hat[i+1] = rk4(y_hat[i], pars, (time[i+1] - time[i]), step_size);
+        Delta_hat[i] = y_hat[i+1] - y_hat[i];
+
+      } else {
+        Delta_hat[i] = DE(y_hat[i], pars)*(time[i] - time[i-1]);
+      }
+    } else {
+      Delta_hat[i] = DE(y_hat[i], pars)*(time[i] - time[i-1]);
+    }
+  }
+}

inst/stan/power_single_ind.stan

@@ -0,0 +1,131 @@
+// Power law model with single individual
+
+functions{
+  //DE function for use with Runge-Kutta method
+  //pars = [ind_coeff, ind_power]
+  real DE(real y, vector pars){
+    return pars[1] * pow((y/pars[3]), -pars[2]);
+  }
+
+  real rk4_step(real y, vector pars, real interval){
+    real k1;
+    real k2;
+    real k3;
+    real k4;
+    real y_hat;
+
+    k1 = DE(y, pars);
+    k2 = DE(y+interval*k1/2.0, pars);
+    k3 = DE(y+interval*k2/2.0, pars);
+    k4 = DE(y+interval*k3, pars);
+
+    y_hat = y + (1.0/6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4) * interval;
+
+    return y_hat;
+  }
+
+  real rk4(real y, vector pars, real interval, real step_size){
+    int steps;
+    real duration;
+    real y_hat;
+    real step_size_temp;
+
+    duration = 0;
+    y_hat = y;
+
+    while(duration < interval){
+      //Determine the relevant step size
+      step_size_temp = min([step_size, interval-duration]);
+
+      //Get next size estimate
+      y_hat = rk4_step(y_hat, pars, step_size_temp);
+
+      //Increment observed duration
+      duration = duration + step_size_temp;
+    }
+
+    return y_hat;
+  }
+}
+
+// Data structure
+data {
+  real step_size;
+  int n_obs;
+  real y_obs[n_obs];
+  int obs_index[n_obs];
+  real time[n_obs];
+  real y_bar;
+  real y_0_obs;
+}
+
+// The parameters accepted by the model.
+parameters {
+  //Individual level
+  real<lower=0> ind_y_0;
+  real<lower=0> ind_coeff;
+  real<lower=0> ind_power;
+
+  //Global level
+  real<lower=0> global_error_sigma;
+}
+
+// The model to be estimated.
+model {
+  real y_hat[n_obs];
+  vector[3] pars;
+
+  pars[1] = ind_coeff;
+  pars[2] = ind_power;
+  pars[3] = y_bar;
+
+  for(i in 1:n_obs){
+
+    if(obs_index[i]==1){//Fits the first size
+      y_hat[i] = ind_y_0;
+    }
+
+    if(i < n_obs){
+      //Estimate next size
+      y_hat[i+1] = rk4(y_hat[i], pars, (time[i+1] - time[i]), step_size);
+    }
+  }
+
+  //Likelihood
+  y_obs ~ normal(y_hat, global_error_sigma);
+
+  //Priors
+  //Individual level
+  ind_y_0 ~ normal(y_0_obs, global_error_sigma);
+  ind_coeff ~lognormal(0, 1);
+  ind_power ~lognormal(0, 1);
+
+  //Global level
+  global_error_sigma ~cauchy(0,5);
+}
+
+generated quantities{
+  real y_hat[n_obs];
+  real Delta_hat[n_obs];
+  vector[3] pars;
+
+  pars[1] = ind_coeff;
+  pars[2] = ind_power;
+  pars[3] = y_bar;
+
+  for(i in 1:n_obs){
+
+    if(obs_index[i]==1){//Fits the first size
+      y_hat[i] = ind_y_0;
+    }
+
+    if(i < n_obs){
+      //Estimate next size
+      y_hat[i+1] = rk4(y_hat[i], pars, (time[i+1] - time[i]), step_size);
+      Delta_hat[i] = y_hat[i+1] - y_hat[i];
+
+    } else {
+      Delta_hat[i] = DE(y_hat[i], pars)*(time[i] - time[i-1]);
+    }
+  }
+}

src/RcppExports.cpp

@@ -17,13 +17,17 @@ RcppExport SEXP _rcpp_module_boot_stan_fit4canham_single_ind_mod();
 RcppExport SEXP _rcpp_module_boot_stan_fit4constant_multi_ind_mod();
 RcppExport SEXP _rcpp_module_boot_stan_fit4constant_single_ind_mod();
 RcppExport SEXP _rcpp_module_boot_stan_fit4linear_mod();
+RcppExport SEXP _rcpp_module_boot_stan_fit4power_multi_ind_mod();
+RcppExport SEXP _rcpp_module_boot_stan_fit4power_single_ind_mod();
 
 static const R_CallMethodDef CallEntries[] = {
     {"_rcpp_module_boot_stan_fit4canham_multi_ind_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4canham_multi_ind_mod, 0},
     {"_rcpp_module_boot_stan_fit4canham_single_ind_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4canham_single_ind_mod, 0},
     {"_rcpp_module_boot_stan_fit4constant_multi_ind_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4constant_multi_ind_mod, 0},
     {"_rcpp_module_boot_stan_fit4constant_single_ind_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4constant_single_ind_mod, 0},
     {"_rcpp_module_boot_stan_fit4linear_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4linear_mod, 0},
+    {"_rcpp_module_boot_stan_fit4power_multi_ind_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4power_multi_ind_mod, 0},
+    {"_rcpp_module_boot_stan_fit4power_single_ind_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4power_single_ind_mod, 0},
     {NULL, NULL, 0}
 };