Constant (#10)

Finalised the testing structure for models with single and multiple individuals, using the constant growth model as the test case. 
- Generating a dataset to be loaded with the model in the fixtures.
- Running the chosen model on the first individual in the simulated data and test that the parameter estimate converges to a reasonable value (with 0.1 of the true value)
- Running the chosen model on multiple individuals and tests whether the output samples are of the right size for the number of chains, iterations, and parameters.

* Notation (#5)

* Updated notation and function structure in the multi individual constant model file for consistency. Updated rmot_run, rmot_models to reflect the changes.

* Removed old constant stan file and cleared out references to it in stanmodels.R, rmot.rmd.

* Updated ignore to exclude compiled binaries.

* Cleared out .o and .so files.

* Added develop branch as a trigger for PR

---------

Co-authored-by: Fonti Kar <[email protected]>

* Added single ind const model, updated rmot_models, rmot_run, rmot vignette. Need to do test files yet.

* Updated the single individual constant model stan file and the testing code.

* Added some comments to the model testing script.

* Expending tests: testing run

* Reorganise if statements

* Rename file to linear

* Added unit testing of model output for linear data.

* Re-added testing structure.

* Found problem with unit testing for const model outputs: object y_single/y_multi in global environment not being seen within rmot_assign_data function.

* Added internal testing data so no data is created within test_that() functions, y_single and y_multi are not evaluated globally #1

* Added internal testing data so no data is created within test_that() functions, y_single and y_multi are not evaluated globally #1

* Skip snapshots on CI

* Updated snaps

* Removed snaps as a testing framework and using expect_equal #1

* Removed skip_on_ci #1

* Updating roxygen version and package doc

* Updated code to generate testing data

* Moved helper code in generating fake data, removed internal testing data

* Updated rmot_assign and tests passing I think

* Completed testing workflow

* Added filepath to saving rds and updated TO DO in assign_data

* Set tolerance

* Updated version for checkout

* Added constant data generation and tests based on model output summary for single and multiple individuals. Tests all complete dynamically and when run in the testing console.

* Init comit for constant branch

* Tess helper funcs

* Regenerate test data using stan seed method

* using rstan::extract instead of summary

* Added constant as trigger

* Moved unique code back to make_constant.R

* Changed test data generation and single- and multi-individual testing structures for constant model.

* Fixed linear regression model testing.

* Resolving issues raised in pull request review: changed the set-up of the constant model and removed reference to the constant branch from R-CMD-checn.yaml.

---------

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

.Rbuildignore

@@ -6,3 +6,6 @@
 ^LICENSE\.md$
 ^\.github$
 ^codecov\.yml$
+^data-raw$
+^tests/testthat/_snaps/$
+

.github/workflows/R-CMD-check.yaml

@@ -4,7 +4,7 @@ on:
   push:
     branches: [main, master]
   pull_request:
-    branches: [main, master]
+    branches: [main, master, develop]
 
 name: R-CMD-check
 
@@ -29,7 +29,7 @@ jobs:
       R_KEEP_PKG_SOURCE: yes
 
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
 
       - uses: r-lib/actions/setup-pandoc@v2
 

.github/workflows/test-coverage.yaml

@@ -4,7 +4,7 @@ on:
   push:
     branches: [main, master]
   pull_request:
-    branches: [main, master]
+    branches: [main, master, develop]
 
 name: test-coverage
 

.gitignore

@@ -51,3 +51,10 @@ rsconnect/
 inst/doc
 /doc/
 /Meta/
+
+# Pre-compiled stan files
+*.o
+*.so
+
+# Bugged directory from snapshot package
+tests/testthat/_snaps

DESCRIPTION

@@ -12,7 +12,7 @@ Description: What the package does (one paragraph).
 License: MIT + file LICENSE
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
 Biarch: true
 Depends: 
     R (>= 3.4.0)
@@ -36,6 +36,7 @@ SystemRequirements: GNU make
 Suggests: 
     knitr,
     rmarkdown,
-    testthat (>= 3.0.0)
+    testthat (>= 3.0.0),
+    withr
 VignetteBuilder: knitr
 Config/testthat/edition: 3

R/rmot-package.R

@@ -15,4 +15,5 @@
 #' @references
 #' Stan Development Team (NA). RStan: the R interface to Stan. R package version 2.26.23. https://mc-stan.org
 #'
+"_PACKAGE"
 NULL

R/rmot_assign_data.R

@@ -10,7 +10,7 @@
 #' rmot_model("linear") |> rmot_assign_data(X = Loblolly$age, Y = Loblolly$height)
 rmot_assign_data <- function(model_template, ...){
   # Grab user expressions
-  user_code <- rlang::enexprs(..., .check_assign = TRUE)
+  user_code <- rlang::enquos(..., .check_assign = TRUE)
 
   # Grab the names
   fields <- names(user_code)
@@ -19,13 +19,14 @@ rmot_assign_data <- function(model_template, ...){
 
   # Evaluate the RHS of expressions (the values)
   data <- purrr::map(user_code,
-                     eval)
+                     ~rlang::eval_tidy(.x, env = rlang::caller_env())
+                     )
 
   for(i in fields){
     model_template <- purrr::list_modify(model_template, !!!data[i])
   }
 
-  #TODO: Check if N is supplied, if not, assign by default to length(X) and give warning
+  #TODO: Check if N is supplied, if not, give error
 
   return(model_template)
 }

R/rmot_models.R

@@ -14,7 +14,8 @@ rmot_model <- function(model=NULL){
 
   output <- switch(model,
                    linear = rmot_lm(),
-                   constant_single = rmot_cgs())
+                   constant_single_ind = rmot_const_single_ind(),
+                   constant_multi_ind = rmot_const_multi_ind())
 
   class(output) <- "rmot_object"
 
@@ -32,19 +33,32 @@ rmot_lm <- function(){
        model = "linear")
 }
 
+#' Data configuration template for constant growth single individual model
+#' @keywords internal
+#' @noRd
+
+rmot_const_single_ind <- function(){
+  list(n_obs = NULL,
+       y_obs = NULL,
+       obs_index = NULL,
+       time = NULL,
+       y_0_obs = NULL,
+       model = "constant_single_ind")
+}
+
 #' Data configuration template for constant growth single species model
 #' @keywords internal
 #' @noRd
 
-rmot_cgs <- function(){
-  list(N_obs = NULL,
-       N_ind = NULL,
-       S_obs = NULL,
-       census = NULL,
-       census_interval = NULL,
-       id_factor = NULL,
-       S_0_obs = NULL,
-       model = "constant_single")
+rmot_const_multi_ind <- function(){
+  list(n_obs = NULL,
+       n_ind = NULL,
+       y_obs = NULL,
+       obs_index = NULL,
+       time = NULL,
+       ind_id = NULL,
+       y_0_obs = NULL,
+       model = "constant_multi_ind")
 }
 
 

R/rmot_run.R

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

R/stanmodels.R

@@ -1,10 +1,11 @@
 # Generated by rstantools.  Do not edit by hand.
 
 # names of stan models
-stanmodels <- c("constant_single", "linear")
+stanmodels <- c("constant_multi_ind", "constant_single_ind", "linear")
 
 # load each stan module
-Rcpp::loadModule("stan_fit4constant_single_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)
 
 # instantiate each stanmodel object

inst/stan/constant_multi_ind.stan

@@ -0,0 +1,96 @@
+//Constant DE - Single species
+functions{
+  //DE function
+  real DE(real beta){
+    return beta;
+  }
+
+  real size_step(real y, real beta, real time){
+    return y + DE(beta) * time;
+  }
+}
+
+// Data structure
+data {
+  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_beta[n_ind];
+
+  real species_beta_mu;
+  real<lower=0> species_beta_sigma;
+
+  //Global level
+  real<lower=0> global_error_sigma;
+}
+
+// The model to be estimated.
+model {
+  real y_hat[n_obs];
+
+  for(i in 1:n_obs){
+    //Fits the first size
+    if(obs_index[i]==1){
+      y_hat[i] = ind_y_0[ind_id[i]];
+    }
+
+    // Estimate next size
+    if(i < n_obs){
+      if(ind_id[i+1]==ind_id[i]){
+        y_hat[i+1] = size_step(y_hat[i], ind_beta[ind_id[i]], (time[i+1]-time[i]));
+      }
+    }
+  }
+
+  //Likelihood
+  y_obs ~ normal(y_hat, global_error_sigma);
+
+  //Priors
+  //Individual level
+  ind_y_0 ~ normal(y_0_obs, global_error_sigma);
+  ind_beta ~ lognormal(species_beta_mu,
+                    species_beta_sigma);
+
+  //Species level
+  species_beta_mu ~ normal(0.1, 1);
+  species_beta_sigma ~cauchy(0.1, 1);
+
+  //Global level
+  global_error_sigma ~cauchy(0.1, 1);
+}
+
+// The output
+generated quantities {
+  real y_hat[n_obs];
+  real Delta_hat[n_obs];
+
+  for(i in 1:n_obs){
+
+    //Fits the first size
+    if(obs_index[i]==1){
+      y_hat[i] = ind_y_0[ind_id[i]];
+    }
+
+    // Estimate next size
+    if(i < n_obs){
+      if(ind_id[i+1]==ind_id[i]){
+        y_hat[i+1] = size_step(y_hat[i], ind_beta[ind_id[i]], (time[i+1]-time[i]));
+        Delta_hat[i] = y_hat[i+1] - y_hat[i];
+      } else {
+        Delta_hat[i] = DE(ind_beta[ind_id[i]]) * (time[i]-time[i-1]);
+      }
+    } else {
+      Delta_hat[i] = DE(ind_beta[ind_id[i]]) * (time[i]-time[i-1]);
+    }
+  }
+}