Merged in master, rm files not needed, work with RMD for README

.gitignore

@@ -31,6 +31,7 @@ vignettes/*.pdf
 # knitr and R markdown default cache directories
 *_cache/
 /cache/
+inst/doc/
 
 # Temporary files created by R markdown
 *.utf8.md
@@ -48,7 +49,6 @@ po/*~
 # RStudio Connect folder
 rsconnect/
 
-inst/doc
 /doc/
 /Meta/
 

DESCRIPTION

@@ -39,7 +39,13 @@ Suggests:
     rmarkdown,
     testthat (>= 3.0.0),
     withr,
-    mnormt
+    mnormt,
+    here,
+    patchwork,
+    deSolve,
+    cowplot,
+    mixtools,
+    MASS
 VignetteBuilder: knitr
 Config/testthat/edition: 3
 LazyData: true

NAMESPACE

@@ -1,15 +1,16 @@
 # Generated by roxygen2: do not edit by hand
 
+export(hmde_affine_de)
 export(hmde_assign_data)
 export(hmde_canham_de)
 export(hmde_const_de)
 export(hmde_extract_estimates)
-export(hmde_linear_de)
 export(hmde_model)
 export(hmde_model_des)
-export(hmde_model_name)
+export(hmde_model_names)
 export(hmde_model_pars)
 export(hmde_plot_de_pieces)
+export(hmde_plot_obs_est_inds)
 export(hmde_run)
 export(hmde_vb_de)
 import(Rcpp)

R/hmde_assign_data.R

@@ -2,21 +2,28 @@
 #'
 #' @param model_template output from hmde_model
 #' @param data Input data tibble with columns including time, y_obs, obs_index, and additionally ind_id for multi-individual models
-#' @param step_size Step size for numerical integration.
 #' @param ... data-masking name-value pairs allowing specific input of elements
 #'
 #' @return updated named list with your data assigned to Stan model parameters
 #' @export
 
-hmde_assign_data <- function(model_template, data = NULL, step_size = NULL, ...){
-  if(!model_template$model %in% hmde_model_name()){
-    stop("Model name not recognised. Run hmde_model_name() to see available models.")
+hmde_assign_data <- function(model_template, data = NULL,...){
+  if(!model_template$model %in% hmde_model_names()){
+    stop("Model name not recognised. Run hmde_model_names() to see available models.")
   }
 
   if(!is.null(data)){ # Use provided tibble
     user_fields <- names(data)
 
+    user_code <- rlang::enquos(..., .check_assign = TRUE)
+    additional_user_fields <- names(user_code)
+    # Evaluate the RHS of expressions (the values)
+    additional_data <- purrr::map(user_code,
+                       ~rlang::eval_tidy(.x, env = rlang::caller_env())
+    )
+
   } else { # Grab user expressions from individual list items and extract data
+    additional_user_fields <- NULL
     user_code <- rlang::enquos(..., .check_assign = TRUE)
     user_fields <- names(user_code)
     # Evaluate the RHS of expressions (the values)
@@ -45,16 +52,21 @@ hmde_assign_data <- function(model_template, data = NULL, step_size = NULL, ...)
   }
 
   for(i in model_fields){ # Iterate through required fields and fill them
-    if(i %in% user_fields){
+    if(i %in% user_fields){ #Check if the user has supplied it in a tibble
       model_template <- purrr::list_modify(model_template, !!!data[i])
-    } else {
+
+    } else if(!is.null(additional_user_fields)){
+      if(i %in% additional_user_fields){ #Check if the user supplied it directly
+        model_template <- purrr::list_modify(model_template, !!!additional_data[i])
+      }
+    }
+
+    if(is.null(model_template[[i]])){ #Catches default tibble transformations
       model_template[[i]] <- switch(
         i,
         n_obs = length(data$y_obs),
         n_ind = length(unique(data$ind_id)),
-        y_0_obs = data$y_obs[which(data$obs_index == 1)],
-        y_bar = mean(data$y_obs),
-        model = model_template$model
+        y_bar = mean(data$y_obs)
       )
     }
 

R/hmde_extract_estimates.R

@@ -1,6 +1,5 @@
 #' Extract samples and return measurement, individual, and population-level estimates
 #'
-#' @param model model name character string
 #' @param fit fitted model Stan fit
 #' @param input_measurement_data data used to fit the model with ind_id, y_obs, time, obs_index tibble
 #'
@@ -9,8 +8,7 @@
 #' @import dplyr
 #' @importFrom stats quantile
 
-hmde_extract_estimates <- function(model = NULL,
-                                   fit = NULL,
+hmde_extract_estimates <- function(fit = NULL,
                                    input_measurement_data = NULL){
   #Check for fit
   if(is.null(fit)){
@@ -21,9 +19,11 @@ hmde_extract_estimates <- function(model = NULL,
     stop("Fit not S4 stanfit type.")
   }
 
+  model <- fit@model_name
   #Check for model
-  if(!model %in% hmde_model_name()){
-    stop("Model name not recognised. Run hmde_model_name() to see available models.")
+  if(!model %in% hmde_model_names()){
+    stop(paste0("Model name not recognised: ", model,
+                " Run hmde_model_names() to see available models."))
   }
 
   #Check for input measurement data
@@ -33,7 +33,7 @@ hmde_extract_estimates <- function(model = NULL,
     }
   }
 
-  estimate_list <- list()
+  estimate_list <- list(model_name = model)
   par_names <- hmde_model_pars(model)
 
   if(grepl("multi", model)){ #Get n_ind for multi-individual

R/hmde_model_des.R

@@ -5,8 +5,8 @@
 #' @export
 
 hmde_model_des <- function(model = NULL){
-  if(!model %in% hmde_model_name()){
-    stop("Model name not recognised. Run hmde_model_name() to see available models.")
+  if(!model %in% hmde_model_names()){
+    stop("Model name not recognised. Run hmde_model_names() to see available models.")
   }
 
   output <- switch(
@@ -17,7 +17,7 @@ hmde_model_des <- function(model = NULL){
     canham_multi_ind = hmde_canham_de,
     vb_single_ind = hmde_vb_de,
     vb_multi_ind = hmde_vb_de,
-    linear_single_ind = hmde_linear_de
+    affine_single_ind = hmde_affine_de
   )
 
   return(output)
@@ -65,15 +65,15 @@ hmde_vb_de <- function(y = NULL, pars = NULL){
     )
 }
 
-#' Differential equation for linear growth single individual model
+#' Differential equation for affine growth single individual model
 #' @param y input real
 #' @param pars list of parameters beta_0, beta_1
 #'
 #' @return value of differential equation at y
 #' @export
 
-hmde_linear_de <- function(y = NULL, pars = NULL){
+hmde_affine_de <- function(y = NULL, pars = NULL){
   return(
-    pars[[1]] + pars[[2]] * y
+    pars[[1]] - pars[[2]] * y
     )
 }

R/hmde_model_names.R

@@ -3,14 +3,14 @@
 #' @return vector of character strings for model names.
 #' @export
 
-hmde_model_name <- function(){
+hmde_model_names <- function(){
   output <- c("constant_single_ind",
               "constant_multi_ind",
               "canham_single_ind",
               "canham_multi_ind",
               "vb_single_ind",
               "vb_multi_ind",
-              "linear_single_ind")
+              "affine_single_ind")
 
   return(output)
 }

R/hmde_model_pars.R

@@ -7,8 +7,8 @@
 
 hmde_model_pars <- function(model=NULL){
 
-  if(!model %in% hmde_model_name()){
-    stop("Model name not recognised. Run hmde_model_name() to see available models.")
+  if(!model %in% hmde_model_names()){
+    stop("Model name not recognised. Run hmde_model_names() to see available models.")
   }
 
   output <- switch(model,
@@ -18,7 +18,7 @@ hmde_model_pars <- function(model=NULL){
                    canham_multi_ind = hmde_canham_multi_ind_pars(),
                    vb_single_ind = hmde_vb_single_ind_pars(),
                    vb_multi_ind = hmde_vb_multi_ind_pars(),
-                   linear_single_ind = hmde_linear_single_ind_pars())
+                   affine_single_ind = hmde_affine_single_ind_pars())
 
   return(output)
 }
@@ -95,13 +95,13 @@ hmde_vb_multi_ind_pars <- function(){
        model = "vb_multi_ind")
 }
 
-#' Parameter names for linear growth single individual model
+#' Parameter names for affine growth single individual model
 #' @keywords internal
 #' @noRd
 #'
-hmde_linear_single_ind_pars <- function(){
+hmde_affine_single_ind_pars <- function(){
   list(measurement_pars_names = c("y_hat"),
        individual_pars_names = c("ind_beta_0", "ind_beta_1"),
-       error_pars_names = c("global_error_sigma"),
-       model = "linear_single_ind")
+       error_pars_names = c(NULL),
+       model = "affine_single_ind")
 }

R/hmde_models.R

@@ -7,8 +7,8 @@
 
 hmde_model <- function(model=NULL){
 
-  if(!model %in% hmde_model_name()){
-    stop("Model name not recognised. Run hmde_model_name() to see available models.")
+  if(!model %in% hmde_model_names()){
+    stop("Model name not recognised. Run hmde_model_names() to see available models.")
   }
 
   output <- switch(model,
@@ -18,7 +18,7 @@ hmde_model <- function(model=NULL){
                    canham_multi_ind = hmde_canham_multi_ind(),
                    vb_single_ind = hmde_vb_single_ind(),
                    vb_multi_ind = hmde_vb_multi_ind(),
-                   linear_single_ind = hmde_linear_single_ind())
+                   affine_single_ind = hmde_affine_single_ind())
 
   class(output) <- "hmde_object"
 
@@ -34,7 +34,8 @@ hmde_const_single_ind <- function(){
        y_obs = NULL,
        obs_index = NULL,
        time = NULL,
-       y_0_obs = NULL,
+       prior_pars_ind_beta = c(0, 2),
+       prior_pars_global_error_sigma = c(0, 2),
        model = "constant_single_ind")
 }
 
@@ -49,7 +50,9 @@ hmde_const_multi_ind <- function(){
        obs_index = NULL,
        time = NULL,
        ind_id = NULL,
-       y_0_obs = NULL,
+       prior_pars_pop_beta_mu = c(0,2),
+       prior_pars_pop_beta_sigma = c(0,2),
+       prior_pars_global_error_sigma = c(0,2),
        model = "constant_multi_ind")
 }
 
@@ -62,7 +65,10 @@ hmde_canham_single_ind <- function(){
        y_obs = NULL,
        obs_index = NULL,
        time = NULL,
-       y_0_obs = NULL,
+       prior_pars_ind_max_growth = c(0,2),
+       prior_pars_ind_size_at_max_growth = c(0,2),
+       prior_pars_ind_k = c(0,2),
+       prior_pars_global_error_sigma = c(0,2),
        model = "canham_single_ind")
 }
 
@@ -77,7 +83,13 @@ hmde_canham_multi_ind <- function(){
        obs_index = NULL,
        time = NULL,
        ind_id = NULL,
-       y_0_obs = NULL,
+       prior_pars_pop_max_growth_mean = c(0,2),
+       prior_pars_pop_max_growth_sd = c(0,2),
+       prior_pars_pop_size_at_max_growth_mean = c(0,2),
+       prior_pars_pop_size_at_max_growth_sd = c(0,2),
+       prior_pars_pop_k_mean = c(0,2),
+       prior_pars_pop_k_sd = c(0,2),
+       prior_pars_global_error_sigma = c(0,2),
        model = "canham_multi_ind")
 }
 
@@ -91,8 +103,10 @@ hmde_vb_single_ind <- function(){
        y_obs = NULL,
        obs_index = NULL,
        time = NULL,
-       y_0_obs = NULL,
        y_bar = NULL,
+       prior_pars_ind_max_size_sd_only = 2,
+       prior_pars_ind_growth_rate = c(0,2),
+       prior_pars_global_error_sigma = c(0,2),
        model = "vb_single_ind")
 }
 
@@ -107,22 +121,28 @@ hmde_vb_multi_ind <- function(){
        obs_index = NULL,
        time = NULL,
        ind_id = NULL,
-       y_0_obs = NULL,
        y_bar = NULL,
+       prior_pars_pop_max_size_mean_sd_only = 2,
+       prior_pars_pop_max_size_sd = c(0,2),
+       prior_pars_pop_growth_rate_mean = c(0,2),
+       prior_pars_pop_growth_rate_sd = c(0,2),
+       prior_pars_global_error_sigma = c(0,2),
        model = "vb_multi_ind")
 }
 
-#' Data configuration template for linear growth single individual model
+#' Data configuration template for affine growth single individual model
 #' @keywords internal
 #' @noRd
 #'
-hmde_linear_single_ind <- function(){
+hmde_affine_single_ind <- function(){
   list(step_size = NULL,
        n_obs = NULL,
        y_obs = NULL,
        obs_index = NULL,
        time = NULL,
-       y_0_obs = NULL,
+       int_method = NULL,
        y_bar = NULL,
-       model = "linear_single_ind")
+       prior_pars_ind_const = c(1,2),
+       prior_pars_ind_beta_1 = c(0,2),
+       model = "affine_single_ind")
 }

R/hmde_plot_de_pieces.R

@@ -27,8 +27,8 @@ hmde_plot_de_pieces <- function(model = NULL,
                                 colour = "#006600",
                                 alpha = 0.4){
   #Check for model
-  if(!model %in% hmde_model_name()){
-    stop("Model name not recognised. Run hmde_model_name() to see available models.")
+  if(!model %in% hmde_model_names()){
+    stop("Model name not recognised. Run hmde_model_names() to see available models.")
   }
 
   if(is.null(model)){
@@ -89,9 +89,9 @@ hmde_ggplot_de_pieces <- function(pars_data,
   plot <- ggplot() +
     xlim(min(y_0), max(y_final)) +
     labs(x = xlab, y = ylab, title = title) +
-    theme_classic() +
-    theme(axis.text=element_text(size=16),
-          axis.title=element_text(size=18,face="bold"))
+    theme_classic()# +
+    #theme(axis.text=element_text(size=16),
+    #      axis.title=element_text(size=18,face="bold"))
 
   for(i in 1:nrow(pars_data)){
     args_list <- list(pars=pars_data[i,-1]) #Remove ind_id