Merge pull request #2 from traitecoevo/list-handling

Constant growth single species template added and an improved way to configure data according for 'model templates'

DESCRIPTION

@@ -15,10 +15,13 @@ Depends:
     R (>= 3.4.0)
 Imports: 
     methods,
+    purrr,
     Rcpp (>= 0.12.0),
     RcppParallel (>= 5.0.1),
+    rlang,
     rstan (>= 2.18.1),
-    rstantools (>= 2.3.1.1)
+    rstantools (>= 2.3.1.1),
+    stats
 LinkingTo: 
     BH (>= 1.66.0),
     Rcpp (>= 0.12.0),

NAMESPACE

@@ -1,6 +1,8 @@
 # Generated by roxygen2: do not edit by hand
 
-export(rmot_lm)
+export(rmot_assign_data)
+export(rmot_model)
+export(rmot_run)
 import(Rcpp)
 import(methods)
 importFrom(RcppParallel,RcppParallelLibs)

R/rmot_assign_data.R

@@ -0,0 +1,31 @@
+#' Assign data to template
+#'
+#' @param model_template output from rmot_model
+#' @param ... data-masking name-value pairs
+#'
+#' @return updated named list with your data assigned to Stan model parameters
+#' @export
+#'
+#' @examples
+#' 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)
+
+  # Grab the names
+  fields <- names(user_code)
+
+  #TODO: Check names are in model_template
+
+  # Evaluate the RHS of expressions (the values)
+  data <- purrr::map(user_code,
+                     eval)
+
+  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
+
+  return(model_template)
+}

R/rmot_models.R

@@ -0,0 +1,51 @@
+#' Select data configuration template for rmot supported model
+#'
+#' @param model model name character string
+#'
+#' @return named list that matches Stan model parameters
+#' @export
+#'
+#' @examples
+#' rmot_model("linear")
+
+rmot_model <- function(model=NULL){
+
+  #TODO: Need a mechanism to check model requested in one that is supported by rmot
+
+  output <- switch(model,
+                   linear = rmot_lm(),
+                   constant_single = rmot_cgs())
+
+  class(output) <- "rmot_object"
+
+  return(output)
+}
+
+#' Data configuration template for linear model
+#' @keywords internal
+#' @noRd
+
+rmot_lm <- function(){
+  list(X = NULL,
+       Y = NULL,
+       N = NULL,
+       model = "linear")
+}
+
+#' 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")
+}
+
+
+

R/rmot_run.R

@@ -0,0 +1,20 @@
+#' Run a linear model in Stan
+#'
+#' @param model_template model template generated by rmot_model and updated by rmot_assign_data
+#' @param ... additional arguments passed to rstan::sampling
+#'
+#' @return Stanfit model output
+#' @export
+#'
+#' @examples
+#' mtcars
+#' rmot_lm(mtcars$mpg, mtcars$disp)
+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, ...))
+
+  return(out)
+}

R/stanmodels.R

@@ -1,10 +1,11 @@
 # Generated by rstantools.  Do not edit by hand.
 
 # names of stan models
-stanmodels <- c("lm")
+stanmodels <- c("constant_single", "linear")
 
 # load each stan module
-Rcpp::loadModule("stan_fit4lm_mod", what = TRUE)
+Rcpp::loadModule("stan_fit4constant_single_mod", what = TRUE)
+Rcpp::loadModule("stan_fit4linear_mod", what = TRUE)
 
 # instantiate each stanmodel object
 stanmodels <- sapply(stanmodels, function(model_name) {

R/utils.R

@@ -0,0 +1,11 @@
+list_rename = function(data, ...) {
+  mapping = sapply(
+    rlang::enquos(...),
+    rlang::as_name
+  )
+  new_names = stats::setNames(nm=names(data))
+  # `new_name = old_name` for consistency with `dplyr::rename`
+  new_names[mapping] = names(mapping)
+  # for `old_name = new_name` use: `new_names[names(mapping)] = mapping`
+  stats::setNames(data, new_names)
+}

inst/stan/constant_single.stan

@@ -0,0 +1,83 @@
+//Constant Growth - Single species
+
+// Data structure
+data {
+  int N_obs;
+  int N_ind;
+  real S_obs[N_obs];
+  int census[N_obs];
+  real census_interval[N_obs];
+  int id_factor[N_obs];
+  real S_0_obs[N_ind];
+}
+
+// The parameters accepted by the model.
+parameters {
+  //Individual level
+  real<lower=0> ind_S_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 S_hat[N_obs];
+  real G_hat[N_obs];
+
+  for(i in 1:N_obs){
+    if(id_factor[i+1]==id_factor[i]){
+      if(census[i]==1){//Fits the first size
+        S_hat[i] = ind_S_0[id_factor[i]];
+      }
+
+      if(i < N_obs){ //Analytic solution
+        G_hat[i] = ind_beta[id_factor[i]];
+        S_hat[i+1] = S_hat[i] + G_hat[i]*census_interval[i+1];
+      }
+    } else {
+      G_hat[i] = 0; //Gives 0 as the growth estimate for the last data point.
+    }
+  }
+
+  //Likelihood
+  S_obs ~ normal(S_hat, global_error_sigma);
+
+  //Priors
+  //Individual level
+  ind_S_0 ~ normal(S_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 S_hat[N_obs];
+  real G_hat[N_obs];
+
+  for(i in 1:N_obs){
+    if(id_factor[i+1]==id_factor[i]){
+      if(census[i]==1){//Fits the first size
+        S_hat[i] = ind_S_0[id_factor[i]];
+      }
+
+      if(i < N_obs){ //Analytic solution
+        G_hat[i] = ind_beta[id_factor[i]];
+        S_hat[i+1] = S_hat[i] + G_hat[i]*census_interval[i+1];
+      }
+    } else {
+      G_hat[i] = 0; //Gives 0 as the growth estimate for the last data point.
+    }
+  }
+}

inst/stan/lm.stan → inst/stan/linear.stan

@@ -1,8 +1,8 @@
 // A linear model
 data {
   int<lower=1> N;
-  vector[N] x;
-  vector[N] y;
+  vector[N] X;
+  vector[N] Y;
 }
 parameters {
   real intercept;
@@ -12,5 +12,5 @@ parameters {
 model {
   // ... priors, etc.
 
-  y ~ normal(intercept + beta * x, sigma);
+  Y ~ normal(intercept + beta * X, sigma);
 }

src/RcppExports.cpp

@@ -12,10 +12,12 @@ Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
 #endif
 
 
-RcppExport SEXP _rcpp_module_boot_stan_fit4lm_mod();
+RcppExport SEXP _rcpp_module_boot_stan_fit4constant_single_mod();
+RcppExport SEXP _rcpp_module_boot_stan_fit4linear_mod();
 
 static const R_CallMethodDef CallEntries[] = {
-    {"_rcpp_module_boot_stan_fit4lm_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4lm_mod, 0},
+    {"_rcpp_module_boot_stan_fit4constant_single_mod", (DL_FUNC) &_rcpp_module_boot_stan_fit4constant_single_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

@@ -0,0 +1,32 @@
+// Generated by rstantools.  Do not edit by hand.
+
+#include <Rcpp.h>
+using namespace Rcpp ;
+#include "stanExports_constant_single.h"
+
+RCPP_MODULE(stan_fit4constant_single_mod) {
+
+
+    class_<rstan::stan_fit<stan_model, boost::random::ecuyer1988> >("rstantools_model_constant_single")
+
+    .constructor<SEXP,SEXP,SEXP>()
+
+
+    .method("call_sampler", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::call_sampler)
+    .method("param_names", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::param_names)
+    .method("param_names_oi", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::param_names_oi)
+    .method("param_fnames_oi", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::param_fnames_oi)
+    .method("param_dims", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::param_dims)
+    .method("param_dims_oi", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::param_dims_oi)
+    .method("update_param_oi", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::update_param_oi)
+    .method("param_oi_tidx", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::param_oi_tidx)
+    .method("grad_log_prob", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::grad_log_prob)
+    .method("log_prob", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::log_prob)
+    .method("unconstrain_pars", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::unconstrain_pars)
+    .method("constrain_pars", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::constrain_pars)
+    .method("num_pars_unconstrained", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::num_pars_unconstrained)
+    .method("unconstrained_param_names", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::unconstrained_param_names)
+    .method("constrained_param_names", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::constrained_param_names)
+    .method("standalone_gqs", &rstan::stan_fit<stan_model, boost::random::ecuyer1988> ::standalone_gqs)
+    ;
+}