Canham model add

DESCRIPTION

@@ -37,6 +37,7 @@ Suggests:
     knitr,
     rmarkdown,
     testthat (>= 3.0.0),
-    withr
+    withr,
+    mnormt
 VignetteBuilder: knitr
 Config/testthat/edition: 3

R/rmot_models.R

@@ -15,7 +15,9 @@ rmot_model <- function(model=NULL){
   output <- switch(model,
                    linear = rmot_lm(),
                    constant_single_ind = rmot_const_single_ind(),
-                   constant_multi_ind = rmot_const_multi_ind())
+                   constant_multi_ind = rmot_const_multi_ind(),
+                   canham_single_ind = rmot_canham_single_ind(),
+                   canham_multi_ind = rmot_canham_multi_ind())
 
   class(output) <- "rmot_object"
 
@@ -61,5 +63,33 @@ rmot_const_multi_ind <- function(){
        model = "constant_multi_ind")
 }
 
+#' Data configuration template for Canham growth single individual model
+#' @keywords internal
+#' @noRd
+
+rmot_canham_single_ind <- function(){
+  list(step_size = NULL,
+       n_obs = NULL,
+       y_obs = NULL,
+       obs_index = NULL,
+       time = NULL,
+       y_0_obs = NULL,
+       model = "canham_single_ind")
+}
+
+#' Data configuration template for Canham growth single species model
+#' @keywords internal
+#' @noRd
 
+rmot_canham_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_0_obs = NULL,
+       model = "canham_multi_ind")
+}
 

R/rmot_run.R

@@ -18,7 +18,9 @@ rmot_run <- function(model_template, ...) {
   out <- switch(model_template$model,
          linear = rstan::sampling(stanmodels$linear, data = 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, ...))
+         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, ...))
 
   return(out)
 }

R/stanmodels.R

@@ -1,9 +1,11 @@
 # Generated by rstantools.  Do not edit by hand.
 
 # names of stan models
-stanmodels <- c("constant_multi_ind", "constant_single_ind", "linear")
+stanmodels <- c("canham_multi_ind", "canham_single_ind", "constant_multi_ind", "constant_single_ind", "linear")
 
 # 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)

inst/stan/canham_multi_ind.stan

@@ -0,0 +1,159 @@
+//Growth function
+functions{
+  //Growth function for use with Runge-Kutta method
+  //pars = (g_max, s_max, k)
+  real DE(real y, vector pars){
+    return pars[1] *
+    exp(-0.5 * pow(log(y / pars[2]) / pars[3], 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_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_max_growth[n_ind];
+  real<lower=0> ind_diameter_at_max_growth[n_ind];
+  real<lower=0> ind_k[n_ind];
+
+  //Species level
+  real species_max_growth_mean;
+  real<lower=0> species_max_growth_sd;
+  real species_diameter_at_max_growth_mean;
+  real<lower=0> species_diameter_at_max_growth_sd;
+  real species_k_mean;
+  real<lower=0> species_k_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_max_growth[ind_id[i]];
+    pars[2] = ind_diameter_at_max_growth[ind_id[i]];
+    pars[3] = ind_k[ind_id[i]];
+
+    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_max_growth ~lognormal(species_max_growth_mean,
+                            species_max_growth_sd);
+  ind_diameter_at_max_growth ~lognormal(species_diameter_at_max_growth_mean,
+                                        species_diameter_at_max_growth_sd);
+  ind_k ~lognormal(species_k_mean,
+                   species_k_sd);
+
+  //Species level
+  species_max_growth_mean ~normal(0, 1);
+  species_max_growth_sd ~cauchy(0, 1);
+  species_diameter_at_max_growth_mean ~normal(0, 1);
+  species_diameter_at_max_growth_sd ~cauchy(0, 1);
+  species_k_mean ~normal(0, 1);
+  species_k_sd ~cauchy(0, 1);
+
+  //Global level
+  global_error_sigma ~cauchy(0, 2);
+}
+
+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_max_growth[ind_id[i]];
+    pars[2] = ind_diameter_at_max_growth[ind_id[i]];
+    pars[3] = ind_k[ind_id[i]];
+
+    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/canham_single_ind.stan

@@ -0,0 +1,132 @@
+//Growth function
+functions{
+  //Growth function for use with Runge-Kutta method
+  //pars = (g_max, s_max, k)
+  real DE(real y, vector pars){
+    return pars[1] *
+    exp(-0.5 * pow(log(y / pars[2]) / pars[3], 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_0_obs;
+}
+
+// The parameters accepted by the model.
+parameters {
+  //Individual level
+  real<lower=0> ind_y_0;
+  real<lower=0> ind_max_growth;
+  real<lower=0> ind_diameter_at_max_growth;
+  real<lower=0> ind_k;
+
+  //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_max_growth;
+  pars[2] = ind_diameter_at_max_growth;
+  pars[3] = ind_k;
+
+  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_max_growth ~lognormal(0, 1);
+  ind_diameter_at_max_growth ~lognormal(3, 1);
+  ind_k ~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_max_growth;
+  pars[2] = ind_diameter_at_max_growth;
+  pars[3] = ind_k;
+
+  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

@@ -12,11 +12,15 @@ Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
 #endif
 
 
+RcppExport SEXP _rcpp_module_boot_stan_fit4canham_multi_ind_mod();
+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();
 
 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},