traitecoevo · Tess-LaCoil · Nov 28, 2024 · Nov 18, 2024 · Nov 18, 2024 · Nov 18, 2024
diff --git a/.Rbuildignore b/.Rbuildignore
@@ -3,3 +3,5 @@
 .github
 LICENSE.md
 codecov.yml
+^doc$
+^Meta$
diff --git a/.gitignore b/.gitignore
@@ -72,3 +72,6 @@ vBtesting.r
 
 # Data prep files
 data-raw/
+
+# Output from vignettes
+figure/
diff --git a/NAMESPACE b/NAMESPACE
@@ -19,5 +19,6 @@ import(methods)
 importFrom(RcppParallel,RcppParallelLibs)
 importFrom(rstan,sampling)
 importFrom(rstantools,rstan_config)
+importFrom(stats,median)
 importFrom(stats,quantile)
 useDynLib(hmde, .registration = TRUE)
diff --git a/R/data.R b/R/data.R
@@ -64,3 +64,19 @@
 #' @references \url{https://doi.org/10.1002/ecy.4140}
 #' @source \url{https://doi.org/10.15146/5xcp-0d46}
 "Tree_Size_Data"
+
+#' Garcinia recondita model estimates - Barro Colorado Island data
+#'
+#' Estimated sizes, individual growth parameters, and population-level hyper-parameters
+#' for Garcinia recondita fit with a Canham growth function hierarchical model.
+#' The data used to fit the model is the Tree_Size_Data object.
+#'
+#' @format ## `Tree_Size_Ests`
+#' A list with 4 elements:
+#' \describe{
+#'   \item{measurement_data}{A tibble with 5 columns that gives information on size observations and estimates.}
+#'   \item{individual_data}{A tibble with 13 columns that gives posterior estimates for individual growth parameters.}
+#'   \item{error_data}{A tibble with 5 columns that gives posterior estimates of the error parameter.}
+#'   \item{population_data}{A tibble with 5 columns that gives posterior estimates for population-level hyper-parameters.}
+#' }
+"Tree_Size_Ests"
diff --git a/R/hmde-package.R b/R/hmde-package.R
@@ -12,6 +12,7 @@
 #' @importFrom rstan sampling
 #' @importFrom rstantools rstan_config
 #' @importFrom RcppParallel RcppParallelLibs
+#' @importFrom stats median
 #'
 #' @references
 #' Stan Development Team (NA). RStan: the R interface to Stan. R package version 2.26.23. https://mc-stan.org

diff --git a/R/hmde_extract_estimates.R b/R/hmde_extract_estimates.R
@@ -33,7 +33,6 @@ hmde_extract_estimates <- function(model = NULL,
     }
   }
 
-
   estimate_list <- list()
   par_names <- hmde_model_pars(model)
 
@@ -72,7 +71,7 @@ hmde_extract_estimates <- function(model = NULL,
   #If model is multi-individual extract population-level estimates and add to list
   if(!is.null(par_names$population_pars_names)){
     estimate_list$population_data <- hmde_extract_pop_par_ests(samples,
-                                                               par_names$population_pars_names)
+                                                               population_pars_names = par_names$population_pars_names)
   }
 
   return(estimate_list)

diff --git a/R/hmde_model_pars.R b/R/hmde_model_pars.R
@@ -28,7 +28,7 @@ hmde_model_pars <- function(model=NULL){
 #' @noRd
 
 hmde_const_single_ind_pars <- function(){
-  list(measurement_pars_names = c("y_hat", "Delta_hat"),
+  list(measurement_pars_names = c("y_hat"),
        individual_pars_names = c("ind_beta"),
        error_pars_names = c("global_error_sigma"),
        model = "constant_single_ind")
@@ -39,7 +39,7 @@ hmde_const_single_ind_pars <- function(){
 #' @noRd
 
 hmde_const_multi_ind_pars <- function(){
-  list(measurement_pars_names = c("y_hat", "Delta_hat"),
+  list(measurement_pars_names = c("y_hat"),
        individual_pars_names = c("ind_beta"),
        population_pars_names = c("pop_beta_mu", "pop_beta_sigma"),
        error_pars_names = c("global_error_sigma"),
@@ -51,8 +51,8 @@ hmde_const_multi_ind_pars <- function(){
 #' @noRd
 
 hmde_canham_single_ind_pars <- function(){
-  list(measurement_pars_names = c("y_hat", "Delta_hat"),
-       individual_pars_names = c("ind_max_growth", "ind_diameter_at_max_growth", "ind_k"),
+  list(measurement_pars_names = c("y_hat"),
+       individual_pars_names = c("ind_max_growth", "ind_size_at_max_growth", "ind_k"),
        error_pars_names = c("global_error_sigma"),
        model = "canham_single_ind")
 }
@@ -62,10 +62,10 @@ hmde_canham_single_ind_pars <- function(){
 #' @noRd
 
 hmde_canham_multi_ind_pars <- function(){
-  list(measurement_pars_names = c("y_hat", "Delta_hat"),
-       individual_pars_names = c("ind_max_growth", "ind_diameter_at_max_growth", "ind_k"),
+  list(measurement_pars_names = c("y_hat"),
+       individual_pars_names = c("ind_max_growth", "ind_size_at_max_growth", "ind_k"),
        population_pars_names = c("pop_max_growth_mean", "pop_max_growth_sd",
-                      "pop_diameter_at_max_growth_mean", "pop_diameter_at_max_growth_sd",
+                      "pop_size_at_max_growth_mean", "pop_size_at_max_growth_sd",
                       "pop_k_mean", "pop_k_sd"),
        error_pars_names = c("global_error_sigma"),
        model = "canham_multi_ind")
@@ -76,7 +76,7 @@ hmde_canham_multi_ind_pars <- function(){
 #' @noRd
 
 hmde_vb_single_ind_pars <- function(){
-  list(measurement_pars_names = c("y_hat", "Delta_hat"),
+  list(measurement_pars_names = c("y_hat"),
        individual_pars_names = c("ind_max_size", "ind_growth_rate"),
        error_pars_names = c("global_error_sigma"),
        model = "vb_single_ind")
@@ -87,7 +87,7 @@ hmde_vb_single_ind_pars <- function(){
 #' @noRd
 
 hmde_vb_multi_ind_pars <- function(){
-  list(measurement_pars_names = c("y_hat", "Delta_hat"),
+  list(measurement_pars_names = c("y_hat"),
        individual_pars_names = c("ind_max_size", "ind_growth_rate"),
        population_pars_names = c("pop_max_size_mean", "pop_max_size_sd",
                       "pop_growth_rate_mean", "pop_growth_rate_sd"),
@@ -100,7 +100,7 @@ hmde_vb_multi_ind_pars <- function(){
 #' @noRd
 #'
 hmde_linear_single_ind_pars <- function(){
-  list(measurement_pars_names = c("y_hat", "Delta_hat"),
+  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")

diff --git a/R/hmde_plot_de_pieces.R b/R/hmde_plot_de_pieces.R
@@ -1,4 +1,6 @@
 #' Plot pieces of chosen differential equation model for each individual.
+#' Structured to take the individual data tibble that is built by the
+#' hmde_extract_estimates function using the ind_par_name_mean estimates.
 #' Function piece will go from the first fitted size to the last.
 #' Accepted ggplot arguments will change the axis labels, title, line colour, alpha
 #'
@@ -23,7 +25,7 @@ hmde_plot_de_pieces <- function(model = NULL,
                                 ylab = "f",
                                 title = NULL,
                                 colour = "#006600",
-                                alpha = 0.2){
+                                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.")
@@ -41,6 +43,13 @@ hmde_plot_de_pieces <- function(model = NULL,
     stop("Measurement data not provided.")
   }
 
+  #Get individual parameter estimates
+  model_par_names <- hmde_model_pars(model)
+  pars_data <- tibble(ind_id = individual_data$ind_id)
+  for(i in model_par_names$individual_pars_names){
+    pars_data[[i]] <- individual_data[[paste0(i, "_mean")]]
+  }
+
   #Extract initial and final sizes for each individual
   initial_and_final_vals <- measurement_data %>%
     group_by(ind_id) %>%
@@ -51,10 +60,10 @@ hmde_plot_de_pieces <- function(model = NULL,
     select(ind_id, y_0, y_final) %>%
     distinct()
 
-  individual_data <- left_join(individual_data, initial_and_final_vals, by="ind_id")
-
   #Generate plot
-  plot <- hmde_ggplot_de_pieces(pars_data = individual_data,
+  plot <- hmde_ggplot_de_pieces(pars_data = pars_data,
+                                y_0 = initial_and_final_vals$y_0,
+                                y_final = initial_and_final_vals$y_final,
                                 DE_function = hmde_model_des(model),
                                 xlab = xlab,
                                 ylab = ylab,
@@ -69,14 +78,16 @@ hmde_plot_de_pieces <- function(model = NULL,
 #' @keywords internal
 #' @noRd
 hmde_ggplot_de_pieces <- function(pars_data,
+                                  y_0,
+                                  y_final,
                                   DE_function,
                                   xlab,
                                   ylab,
                                   title,
                                   colour,
                                   alpha){
   plot <- ggplot() +
-    xlim(min(pars_data$y_0), max(pars_data$y_final)) +
+    xlim(min(y_0), max(y_final)) +
     labs(x = xlab, y = ylab, title = title) +
     theme_classic() +
     theme(axis.text=element_text(size=16),
@@ -86,8 +97,8 @@ hmde_ggplot_de_pieces <- function(pars_data,
     args_list <- list(pars=pars_data[i,-1]) #Remove ind_id
     plot <- plot +
       geom_function(fun=DE_function, args=args_list,
-                    colour=colour, linewidth=1,
-                    xlim=c(pars_data$y_0[i], pars_data$y_final[i]))
+                    colour=colour, linewidth=1, alpha = alpha,
+                    xlim=c(y_0[i], y_final[i]))
   }
 
   return(plot)

diff --git a/data/Lizard_Size_Data.rda b/data/Lizard_Size_Data.rda
diff --git a/data/Tree_Size_Ests.rda b/data/Tree_Size_Ests.rda
diff --git a/inst/stan/vb_multi_ind.stan b/inst/stan/vb_multi_ind.stan
@@ -99,7 +99,7 @@ generated quantities{
         //Estimate next size
         y_hat[i+1] = solution(time[i+1], pars) + y_bar;
        Delta_hat[i] = y_hat[i+1] - y_hat[i];
-      } else { #Estimate next growth based on same time to last.
+      } else {// Estimate next growth based on same time to last.
         temp_y_final = solution(2*time[i] - time[i-1], pars) + y_bar;
         Delta_hat[i] = temp_y_final - y_hat[i];
       }

diff --git a/inst/stan/vb_single_ind.stan b/inst/stan/vb_single_ind.stan
@@ -84,7 +84,7 @@ generated quantities{
       y_hat[i+1] = solution(time[i+1], pars) + y_bar;
       Delta_hat[i] = y_hat[i+1] - y_hat[i];
 
-    } else { #Estimate next growth based on same time to last.
+    } else { // Estimate next growth based on same time to last.
       temp_y_final = solution(2*time[i] - time[i-1], pars) + y_bar;
       Delta_hat[i] = temp_y_final - y_hat[i];
     }

diff --git a/man/Tree_Size_Ests.Rd b/man/Tree_Size_Ests.Rd
diff --git a/man/hmde_plot_de_pieces.Rd b/man/hmde_plot_de_pieces.Rd
diff --git a/tests/testthat/test-hmde_inbuilt_data.R b/tests/testthat/test-hmde_inbuilt_data.R
@@ -3,16 +3,19 @@ test_that("Datasets: existence", {
   expect_visible(Lizard_Size_Data)
   expect_visible(Trout_Size_Data)
   expect_visible(Tree_Size_Data)
+  expect_visible(Tree_Size_Ests)
 })
 
 test_that("Datasets: size", {
   expect_equal(object = c(nrow(Lizard_Size_Data),
                           ncol(Lizard_Size_Data)),
-               expected = c(332, 4))
+               expected = c(328, 4))
   expect_equal(object = c(nrow(Trout_Size_Data),
                           ncol(Trout_Size_Data)),
                expected = c(135, 4))
   expect_equal(object = c(nrow(Tree_Size_Data),
                           ncol(Tree_Size_Data)),
                expected = c(300, 4))
+  expect_equal(object = length(Tree_Size_Ests),
+               expected = 4)
 })
diff --git a/vignettes/.gitignore b/vignettes/.gitignore
@@ -1,2 +1,3 @@
 *.html
 *.R
+*.Rhistory