Constant

.Rbuildignore

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

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

@@ -2,9 +2,9 @@
 # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
 on:
   push:
-    branches: [main, master]
+    branches: [main, master, constant]
   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,97 @@
+//Constant Growth - Single species
+functions{
+  //Growth function
+  real growth(real beta){
+    return beta;
+  }
+}
+
+// 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];
+  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 growth rate
+    Delta_hat[i] = growth(ind_beta[ind_id[i]]);
+
+    // Update next size if relevant
+    if(i < n_obs){
+      if(ind_id[i+1]==ind_id[i]){
+        //Analytic solution
+        y_hat[i+1] = y_hat[i] + Delta_hat[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 growth rate
+    Delta_hat[i] = growth(ind_beta[ind_id[i]]);
+
+    // Update next size if relevant
+    if(i < n_obs){
+      if(ind_id[i+1]==ind_id[i]){
+        //Analytic solution
+        y_hat[i+1] = y_hat[i] + Delta_hat[i]*(time[i+1]-time[i]);
+      }
+    }
+  }
+}