Added a function to plot sizes over time in parallel with the DE plot…

NAMESPACE

@@ -10,6 +10,7 @@ export(hmde_model_des)
 export(hmde_model_name)
 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_plot_obs_est_inds.R

@@ -0,0 +1,78 @@
+#' Plot estimated and observed values over time for a chosen number of individuals based
+#' on posterior estimates. Structured to take in the measurement_data tibble constructed by
+#' the hmde_extract_estimates function.
+#'
+#' @param measurement_data tibble with estimated measurements
+#' @param ind_id_vec vector with list of ind_id values
+#' @param n_ind_to_plot integer giving number of individuals to plot if not speciried
+#' @param xlab character string for replacement x axis label
+#' @param ylab character string for replacement y axis label
+#' @param title character string for replacement plot title
+#'
+#' @return ggplot object
+#' @export
+#' @import ggplot2
+#' @import dplyr
+
+hmde_plot_obs_est_inds <- function(ind_id_vec = NULL,
+                                   n_ind_to_plot = NULL,
+                                   measurement_data = NULL,
+                                   xlab = "Time",
+                                   ylab = "Y(t)",
+                                   title = NULL){
+  if(is.null(measurement_data)){
+    stop("Measurement data not provided.")
+  }
+
+  if(!is.null(ind_id_vec)){
+    for(i in ind_id_vec){ #Error if ID nujmber not in data.
+      if(!i %in% measurement_data$ind_id){
+        stop(paste0("Ind ID values not recognised: ", i))
+      }
+    }
+
+    plot_data <- measurement_data %>%
+      filter(ind_id %in% ind_id_vec)
+
+  } else if(is.null(n_ind_to_plot)){ #Error if no info for which inds to plot
+    stop("Neither ind. ID values nor a number of individuals provided.")
+
+  } else if(length(unique(measurement_data$ind_id)) < n_ind_to_plot){
+    stop("Number of individuals to plot larger than sample size, please provide a smaller number.")
+
+  } else {
+    sample_ids <- sample(unique(measurement_data$ind_id), size=n_ind_to_plot)
+    plot_data <- measurement_data %>%
+      filter(ind_id %in% sample_ids)
+  }
+
+  #Generate plot
+  plot <- hmde_ggplot_obs_est_inds(plot_data = plot_data,
+                                xlab = xlab,
+                                ylab = ylab,
+                                title = title)
+
+  return(plot)
+}
+
+#' Produce plot for plot_obs_est_inds
+#' @keywords internal
+#' @noRd
+hmde_ggplot_obs_est_inds <- function(plot_data,
+                                     xlab,
+                                     ylab,
+                                     title){
+  plot <- ggplot(data=plot_data, aes(group = ind_id)) +
+    geom_point(aes(x = time, y=y_obs, colour = as.factor(ind_id)),
+               shape = 1) +
+    geom_line(aes(x = time, y=y_obs, colour = as.factor(ind_id)),
+              linetype = "dashed") +
+    geom_point(aes(x = time, y=y_hat, colour = as.factor(ind_id)),
+               shape = 2) +
+    geom_line(aes(x = time, y=y_hat, colour = as.factor(ind_id)),
+              linetype = "solid") +
+    labs(x=xlab, y=ylab, colour="Ind. ID") +
+    theme_classic()
+
+  return(plot)
+}

README.md

@@ -12,7 +12,8 @@ coverage](https://codecov.io/gh/traitecoevo/hmde/branch/master/graph/badge.svg)]
 experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental)
 <!-- badges: end -->
 
-The goal of hmde is to implement hierarchical Bayesian longitudinal models to solve the Bayesian inverse problem of estimating differential equation parameters based on repeat measurement surveys. Estimation is done using Markov Chain Monte Carlo, implemented through
+The goal of hmde is to fit a model for the rate of change in some quantity based on a set of pre-defined functions arising from ecological applications. We estimate differential equation parameters from repeated observations of a process, such as growth rate parameters to data of sizes over time.
+In other language, `hmde` implements hierarchical Bayesian longitudinal models to solve the Bayesian inverse problem of estimating differential equation parameters based on repeat measurement surveys. Estimation is done using Markov Chain Monte Carlo, implemented through
 [Stan](https://mc-stan.org/) via [RStan](https://mc-stan.org/users/interfaces/rstan), built under [R](https://cran.r-project.org/) 4.3.3. The inbuilt models are based on case studies in ecology.
 
 ## The Maths

tests/testthat/test-hmde_plot_obs_est_inds.R

@@ -0,0 +1,30 @@
+test_that("Execution and output: plot_obs_est_inds function", {
+  plot <- hmde_plot_obs_est_inds(n_ind_to_plot = 2,
+                                 measurement_data = Tree_Size_Ests$measurement_data)
+
+  expect_named(plot)
+
+  expect_visible(plot)
+
+  expect_type(plot, "list")
+})
+
+
+test_that("Execution and output: bad input", {
+  expect_error(
+    hmde_plot_obs_est_inds(measurement_data = Tree_Size_Ests$measurement_data)
+  )
+
+  expect_error(
+    hmde_plot_obs_est_inds()
+  )
+
+  expect_error(
+    hmde_plot_obs_est_inds(individual_data = Tree_Size_Ests$measurement_data)
+  )
+
+  expect_error(
+    hmde_plot_obs_est_inds(n_ind_to_plot = 10^3,
+      measurement_data = Tree_Size_Ests$measurement_data)
+  )
+})

vignettes/canham.Rmd

@@ -120,21 +120,8 @@ hist(measurement_data_transformed$est_growth_rate,
 cor(measurement_data_transformed$y_obs, measurement_data_transformed$y_hat)^2
 
 #Plots of size over time for a sample of 5 individuals
-sample_ids <- sample(1:nrow(Tree_Size_Ests$individual_data), size=5)
-plot_data <- measurement_data_transformed %>%
-  filter(ind_id %in% sample_ids)
-
-ggplot(data=plot_data, aes(group = ind_id)) +
-  geom_point(aes(x = time, y=y_obs, colour = as.factor(ind_id)),
-             shape = 1) +
-  geom_line(aes(x = time, y=y_obs, colour = as.factor(ind_id)),
-            linetype = "dashed") +
-  geom_point(aes(x = time, y=y_hat, colour = as.factor(ind_id)),
-             shape = 2) +
-  geom_line(aes(x = time, y=y_hat, colour = as.factor(ind_id)),
-            linetype = "solid") +
-  labs(x="Time (years)", y="DBH (cm)", colour="Ind. ID") +
-  theme_classic()
+hmde_plot_obs_est_inds(n_ind_to_plot = 5,
+                       measurement_data = Tree_Size_Ests$measurement_data)
 ```
 
 Individual parameter analysis, growth function plots follow this.

vignettes/constant-growth.Rmd

@@ -300,22 +300,8 @@ In the next block, we are looking at 5 random individuals to start with because
 cor(measurement_data_transformed$y_obs, measurement_data_transformed$y_hat)^2
 
 #Plots of size over time for a sample of 5 individuals
-sample_ids <- sample(1:nrow(trout_constant_estimates$individual_data), size=5)
-
-# Line plots of sizes over time
-measurement_data_transformed %>%
-  filter(ind_id %in% sample_ids) %>%
-  ggplot(aes(group = ind_id)) +
-  geom_point(aes(x = time, y=y_obs, colour = as.factor(ind_id)), 
-             shape = 1) +
-  geom_line(aes(x = time, y=y_obs, colour = as.factor(ind_id)), 
-            linetype = "dashed") +
-  geom_point(aes(x = time, y=y_hat, colour = as.factor(ind_id)), 
-             shape = 2) +
-  geom_line(aes(x = time, y=y_hat, colour = as.factor(ind_id)), 
-            linetype = "solid") +
-  labs(x="Time (years)", y="Size (cm)", colour="Ind. ID") +
-  theme_classic()
+hmde_plot_obs_est_inds(n_ind_to_plot = 5,
+                       measurement_data = trout_constant_estimates$measurement_data)
 ```
 
 ### Indvidual growth functions ($\beta$)

vignettes/von-bertalanffy.Rmd

@@ -133,21 +133,8 @@ hist(measurement_data_transformed$est_growth_rate,
 cor(measurement_data_transformed$y_obs, measurement_data_transformed$y_hat)^2
 
 #Plots of size over time for a sample of 5 individuals
-sample_ids <- sample(1:nrow(lizard_vb_estimates$individual_data), size=5)
-plot_data <- measurement_data_transformed %>%
-  filter(ind_id %in% sample_ids)
-
-ggplot(data=plot_data, aes(group = ind_id)) +
-  geom_point(aes(x = time, y=y_obs, colour = as.factor(ind_id)), 
-             shape = 1) +
-  geom_line(aes(x = time, y=y_obs, colour = as.factor(ind_id)), 
-            linetype = "dashed") +
-  geom_point(aes(x = time, y=y_hat, colour = as.factor(ind_id)), 
-             shape = 2) +
-  geom_line(aes(x = time, y=y_hat, colour = as.factor(ind_id)), 
-            linetype = "solid") +
-  labs(x="Time (days)", y="Size (mm)", colour="Ind. ID") +
-  theme_classic()
+hmde_plot_obs_est_inds(n_ind_to_plot = 5,
+                       measurement_data = lizard_vb_estimates$measurement_data)
 ```
 
 We have two parameters at the individual level and are interested in both their separate distributions, and if we see evidence of a relationship between them. We can also use the individual parameter estimates and estimated sizes to plot the growth function pieces.