Integration update (#40)

* First attempt at implementing analytic solution for vB model single individual.

* Corrected compile errors.

* Added analytic solution for log-transformed power law model.

* Updated linear/log-linear models.

* Updated parameter names and finalised vB model.

* Updated power model to correct errors in y_bar transformation and time indexing. Updated data generation to have optional step size argument for numerical projection.

* Updated vignettes to get started on the one for mathematicians.

* Removed draft vignette to put in their own branch.

* Added separate stan files for a numeric and analytic solution as something is going wrong in the analytic solution. Need to check the maths.

* Updated DE plotting function to better handle parameter estimates.

* Updated power model.

* Updated lizard and tree data, the former to be length rather than mass, the latter to be G. recondita.

* Updated integratio for Canham to remove the RK4 solver and use the stiff ODE solver from stan.

* Swapped to RK45 to test for stiffness and compilation errors.

* Took out step_size argument from canham.

* Updated testing.

* Updates to testing again.

R/data.R

@@ -2,6 +2,7 @@
 #'
 #' A subset of data from the SUSTAIN trout capture-recapture data set
 #' from Moe et al. (2020).
+#' Observations are of total body length in centimetres.
 #' Data prepared by taking a stratified sample of individual IDs based on the
 #' number of observations per individual: 25 individuals with 2 observations,
 #' 15 with 3, 10 with 4. Within the groups a simple random sample
@@ -24,21 +25,22 @@
 #'
 #' A subset of data from Kar, Nakagawa, and Noble (2024), used to model growth
 #' behaviour in a skink species.
+#' Observations are of the length from the tip of the nose to the start of the cloaca.
 #' Data was prepared by taking a simple random sample with replacement of 50
 #' individual IDs among individuals with at least 5 observations each. Data was
 #' then transformed to conform to the needs of a model data set in the package.
 #'
-#' @format ## `Lizard_Mass_Data`
+#' @format ## `Lizard_Size_Data`
 #' A data frame with 336 rows and 4 columns:
 #' \describe{
 #'   \item{ind_id}{ID number for individual}
 #'   \item{time}{Days since first observation.}
-#'   \item{y_obs}{Individual mass in kilograms.}
+#'   \item{y_obs}{Individual size in mm.}
 #'   \item{obs_index}{Index of observations for individual}
 #' }
 #'
 #' @source \url{https://osf.io/hjkxd/}
-"Lizard_Mass_Data"
+"Lizard_Size_Data"
 
 #' Protium tenufolium - Barro Colorado Island data
 #'
@@ -63,14 +65,14 @@
 #' @source \url{https://doi.org/10.15146/5xcp-0d46}
 "Small_Tree_Size_Data"
 
-#' Hirtella triandra - Barro Colorado Island data
+#' Garcinia recondita - Barro Colorado Island data
 #'
 #' A subset of data from the Barro Colorado Island long term forest plot
 #' managed by the Smithsonian Tropical Research Institute (Condit et al. 2019).
 #' Data was prepared by taking a simple random sample without replacement of
-#' 30 individual IDs from Hirtella triandra. The sampling frame was restricted
+#' 30 individual IDs from Garcinia recondita. The sampling frame was restricted
 #' to individuals with 6 observations since 1990, and a difference between
-#' observed max and min sizes of more than 7cm in order to avoid identifiability
+#' observed first and last sizes of more than 3cm in order to avoid identifiability
 #' issues. Data was then transformed and renamed to match the required structure
 #' to act as demonstration for the package.
 #'

R/hmde_assign_data.R

@@ -50,7 +50,6 @@ hmde_assign_data <- function(model_template, data = NULL, step_size = NULL, ...)
     } else {
       model_template[[i]] <- switch(
         i,
-        step_size = step_size,
         n_obs = length(data$y_obs),
         n_ind = length(unique(data$ind_id)),
         y_0_obs = data$y_obs[which(data$obs_index == 1)],

R/hmde_model_des.R

@@ -70,7 +70,7 @@ hmde_power_de <- function(y = NULL, pars = NULL){
 
 #' Differential equation for von Bertalanffy growth single and multi- individual models
 #' @param y input real
-#' @param pars list of parameters beta, Y_max
+#' @param pars list of parameteters Y_max, growth_rate
 #'
 #' @return value of differential equation at y
 #' @export

R/hmde_model_pars.R

@@ -66,9 +66,9 @@ hmde_canham_single_ind_pars <- function(){
 hmde_canham_multi_ind_pars <- function(){
   list(measurement_pars_names = c("y_hat", "Delta_hat"),
        individual_pars_names = c("ind_max_growth", "ind_diameter_at_max_growth", "ind_k"),
-       population_pars_names = c("pop_max_growth_mu", "pop_max_growth_sigma",
-                      "pop_diameter_at_max_growth_mu", "pop_diameter_at_max_growth_sigma",
-                      "pop_k_mu", "pop_k_sigma"),
+       population_pars_names = c("pop_max_growth_mean", "pop_max_growth_sd",
+                      "pop_diameter_at_max_growth_mean", "pop_diameter_at_max_growth_sd",
+                      "pop_k_mean", "pop_k_sd"),
        error_pars_names = c("global_error_sigma"),
        model = "canham_multi_ind")
 }
@@ -91,8 +91,8 @@ hmde_power_single_ind_pars <- function(){
 hmde_power_multi_ind_pars <- function(){
   list(measurement_pars_names = c("y_hat", "Delta_hat"),
        individual_pars_names = c("ind_coeff", "ind_power"),
-       population_pars_names = c("pop_coeff_mu", "pop_coeff_sigma",
-                      "pop_power_mu", "pop_power_sigma"),
+       population_pars_names = c("pop_coeff_mean", "pop_coeff_sd",
+                      "pop_power_mean", "pop_power_sd"),
        error_pars_names = c("global_error_sigma"),
        model = "power_multi_ind")
 }
@@ -115,8 +115,8 @@ hmde_vb_single_ind_pars <- function(){
 hmde_vb_multi_ind_pars <- function(){
   list(measurement_pars_names = c("y_hat", "Delta_hat"),
        individual_pars_names = c("ind_max_size", "ind_growth_rate"),
-       population_pars_names = c("pop_max_size_mu", "pop_max_size_sigma",
-                      "pop_growth_rate_mu", "pop_growth_rate_sigma"),
+       population_pars_names = c("pop_max_size_mean", "pop_max_size_sd",
+                      "pop_growth_rate_mean", "pop_growth_rate_sd"),
        error_pars_names = c("global_error_sigma"),
        model = "vb_multi_ind")
 }

R/hmde_models.R

@@ -60,8 +60,7 @@ hmde_const_multi_ind <- function(){
 #' @noRd
 
 hmde_canham_single_ind <- function(){
-  list(step_size = NULL,
-       n_obs = NULL,
+  list(n_obs = NULL,
        y_obs = NULL,
        obs_index = NULL,
        time = NULL,
@@ -74,8 +73,7 @@ hmde_canham_single_ind <- function(){
 #' @noRd
 
 hmde_canham_multi_ind <- function(){
-  list(step_size = NULL,
-       n_obs = NULL,
+  list(n_obs = NULL,
        n_ind = NULL,
        y_obs = NULL,
        obs_index = NULL,
@@ -90,8 +88,7 @@ hmde_canham_multi_ind <- function(){
 #' @noRd
 
 hmde_power_single_ind <- function(){
-  list(step_size = NULL,
-       n_obs = NULL,
+  list(n_obs = NULL,
        y_obs = NULL,
        obs_index = NULL,
        time = NULL,
@@ -105,8 +102,7 @@ hmde_power_single_ind <- function(){
 #' @noRd
 
 hmde_power_multi_ind <- function(){
-  list(step_size = NULL,
-       n_obs = NULL,
+  list(n_obs = NULL,
        n_ind = NULL,
        y_obs = NULL,
        obs_index = NULL,
@@ -122,8 +118,7 @@ hmde_power_multi_ind <- function(){
 #' @noRd
 
 hmde_vb_single_ind <- function(){
-  list(step_size = NULL,
-       n_obs = NULL,
+  list(n_obs = NULL,
        y_obs = NULL,
        obs_index = NULL,
        time = NULL,
@@ -137,8 +132,7 @@ hmde_vb_single_ind <- function(){
 #' @noRd
 
 hmde_vb_multi_ind <- function(){
-  list(step_size = NULL,
-       n_obs = NULL,
+  list(n_obs = NULL,
        n_ind = NULL,
        y_obs = NULL,
        obs_index = NULL,

R/hmde_plot_de_pieces.R

@@ -83,7 +83,7 @@ hmde_ggplot_de_pieces <- function(pars_data,
           axis.title=element_text(size=18,face="bold"))
 
   for(i in 1:nrow(pars_data)){
-    args_list <- list(pars=pars_data[i,])
+    args_list <- list(pars=pars_data[i,-1]) #Remove ind_id
     plot <- plot +
       geom_function(fun=DE_function, args=args_list,
                     colour=colour, linewidth=1,

R/stanmodels.R

@@ -1,7 +1,7 @@
 # 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_single_ind", "power_multi_ind", "power_single_ind", "vb_multi_ind", "vb_single_ind")
+stanmodels <- c("canham_multi_ind", "canham_single_ind", "constant_multi_ind", "constant_single_ind", "linear_single_ind", "power_multi_ind", "power_single_ind", "power_single_ind_analytic", "power_single_ind_numeric", "vb_multi_ind", "vb_single_ind")
 
 # load each stan module
 Rcpp::loadModule("stan_fit4canham_multi_ind_mod", what = TRUE)
@@ -11,6 +11,8 @@ Rcpp::loadModule("stan_fit4constant_single_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4linear_single_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4power_multi_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4power_single_ind_mod", what = TRUE)
+Rcpp::loadModule("stan_fit4power_single_ind_analytic_mod", what = TRUE)
+Rcpp::loadModule("stan_fit4power_single_ind_numeric_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4vb_multi_ind_mod", what = TRUE)
 Rcpp::loadModule("stan_fit4vb_single_ind_mod", what = TRUE)
 

inst/stan/canham_multi_ind.stan

@@ -1,56 +1,15 @@
 //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;
+  vector DE(real t, vector y, real max_growth, real size_at_max, real k){
+    vector[size(y)] dydt = max_growth *
+    exp(-0.5 * square(log(y / size_at_max) / k));
+    return dydt;
   }
 }
 
 // Data structure
 data {
-  real step_size;
   int n_obs;
   int n_ind;
   real y_obs[n_obs];
@@ -83,22 +42,22 @@ parameters {
 // The model to be estimated.
 model {
   real y_hat[n_obs];
-  vector[3] pars;
+  vector[1] y_temp;
 
   for(i in 1:n_obs){
-    // Initialise the parameters for the observation
-    pars[1] = ind_max_growth[ind_id[i]];
-    pars[2] = ind_size_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]){
+        y_temp[1] = y_hat[i];
         //Estimate next size
-        y_hat[i+1] = rk4(y_hat[i], pars, (time[i+1] - time[i]), step_size);
+        y_hat[i+1] = ode_rk45(DE, y_temp,
+          time[i], {time[i+1]},
+          ind_max_growth[ind_id[i]],
+          ind_size_at_max_growth[ind_id[i]],
+          ind_k[ind_id[i]])[1][1];
       }
     }
   }
@@ -130,30 +89,24 @@ model {
 
 generated quantities{
   real y_hat[n_obs];
-  real Delta_hat[n_obs];
-  vector[3] pars;
+  vector[1] y_temp;
 
   for(i in 1:n_obs){
-    // Initialise the parameters for the observation
-    pars[1] = ind_max_growth[ind_id[i]];
-    pars[2] = ind_size_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]){
+        y_temp[1] = y_hat[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]);
+        y_hat[i+1] = ode_rk45(DE, y_temp,
+          time[i], {time[i+1]},
+          ind_max_growth[ind_id[i]],
+          ind_size_at_max_growth[ind_id[i]],
+          ind_k[ind_id[i]])[1][1];
       }
-    } else {
-      Delta_hat[i] = DE(y_hat[i], pars)*(time[i] - time[i-1]);
     }
   }
 }