Constant single individual

.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
 

.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,7 @@ rsconnect/
 inst/doc
 /doc/
 /Meta/
+
+# Pre-compiled stan files
+*.o
+*.so

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,89 @@
+//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){
+    if(ind_id[i+1]==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){ //Analytic solution
+        Delta_hat[i] = growth(ind_beta[ind_id[i]]);
+        y_hat[i+1] = y_hat[i] + Delta_hat[i]*(time[i+1]-time[i]);
+      }
+    } else {
+      Delta_hat[i] = 0; //Gives 0 as the growth estimate for the last data point.
+    }
+  }
+
+  //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){
+    if(ind_id[i+1]==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){ //Analytic solution
+        Delta_hat[i] = growth(ind_beta[ind_id[i]]);
+        y_hat[i+1] = y_hat[i] + Delta_hat[i]*(time[i+1]-time[i]);
+      }
+    } else {
+      Delta_hat[i] = 0; //Gives 0 as the growth estimate for the last data point.
+    }
+  }
+}

inst/stan/constant_single_ind.stan

@@ -0,0 +1,75 @@
+//Constant Growth - Single species
+functions{
+  //Growth function
+  real growth(real y, real beta){
+    return beta;
+  }
+}
+
+// Data structure
+data {
+  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_beta;
+
+  //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){
+    if(obs_index[i]==1){//Fits the first size
+      y_hat[i] = ind_y_0;
+    }
+
+    if(i < n_obs){ //Analytic solution
+      Delta_hat[i] = growth(y_hat[i], ind_beta);
+      y_hat[i+1] = y_hat[i] + Delta_hat[i]*(time[i+1]-time[i]);
+    } else {
+      Delta_hat[i] = 0;
+    }
+  }
+
+  //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(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){
+    if(obs_index[i]==1){//Fits the first size
+      y_hat[i] = ind_y_0;
+    }
+
+    if(i < n_obs){ //Analytic solution
+      Delta_hat[i] = growth(y_hat[i], ind_beta);
+      y_hat[i+1] = y_hat[i] + Delta_hat[i]*(time[i+1]-time[i]);
+    } else {
+      Delta_hat[i] = 0;
+    }
+  }
+}

rmot.Rproj

@@ -19,4 +19,3 @@ LineEndingConversion: Posix
 BuildType: Package
 PackageUseDevtools: Yes
 PackageInstallArgs: --no-multiarch --with-keep.source
-PackageRoxygenize: rd,collate,namespace

src/RcppExports.cpp

@@ -12,11 +12,13 @@ Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
 #endif
 
 
-RcppExport SEXP _rcpp_module_boot_stan_fit4constant_single_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_fit4constant_single_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4constant_single_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},
     {NULL, NULL, 0}
 };

src/stanExports_constant_single.cc → src/stanExports_constant_multi_ind.cc

@@ -2,12 +2,12 @@
 
 #include <Rcpp.h>
 using namespace Rcpp ;
-#include "stanExports_constant_single.h"
+#include "stanExports_constant_multi_ind.h"
 
-RCPP_MODULE(stan_fit4constant_single_mod) {
+RCPP_MODULE(stan_fit4constant_multi_ind_mod) {
 
 
-    class_<rstan::stan_fit<stan_model, boost::random::ecuyer1988> >("rstantools_model_constant_single")
+    class_<rstan::stan_fit<stan_model, boost::random::ecuyer1988> >("rstantools_model_constant_multi_ind")
 
     .constructor<SEXP,SEXP,SEXP>()