From 466f0887609fbacba40440c6179e6cf4a85911fa Mon Sep 17 00:00:00 2001 From: Tess-LaCoil Date: Wed, 4 Dec 2024 15:28:10 +1100 Subject: [PATCH 1/3] Added a function to plot sizes over time in parallel with the DE plot function. Changed wording in README to be more accessible. --- R/hmde_plot_obs_est_inds.R | 81 ++++++++++++++++++++ README.md | 3 +- tests/testthat/test-hmde_plot_obs_est_inds.R | 30 ++++++++ vignettes/canham.Rmd | 17 +--- vignettes/constant-growth.Rmd | 18 +---- vignettes/von-bertalanffy.Rmd | 17 +--- 6 files changed, 119 insertions(+), 47 deletions(-) create mode 100644 R/hmde_plot_obs_est_inds.R create mode 100644 tests/testthat/test-hmde_plot_obs_est_inds.R diff --git a/R/hmde_plot_obs_est_inds.R b/R/hmde_plot_obs_est_inds.R new file mode 100644 index 0000000..1e381a4 --- /dev/null +++ b/R/hmde_plot_obs_est_inds.R @@ -0,0 +1,81 @@ +#' 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 model model name character string +#' @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 +#' @param colour character string for replacement line colour +#' @param alpha real number for replacement alpha value +#' +#' @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) +} diff --git a/README.md b/README.md index ceeaa73..b7c7356 100644 --- a/README.md +++ b/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) -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 diff --git a/tests/testthat/test-hmde_plot_obs_est_inds.R b/tests/testthat/test-hmde_plot_obs_est_inds.R new file mode 100644 index 0000000..c98b1af --- /dev/null +++ b/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) + ) +}) diff --git a/vignettes/canham.Rmd b/vignettes/canham.Rmd index fe32bd0..8a051e2 100644 --- a/vignettes/canham.Rmd +++ b/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. diff --git a/vignettes/constant-growth.Rmd b/vignettes/constant-growth.Rmd index 74d8669..b62733a 100644 --- a/vignettes/constant-growth.Rmd +++ b/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$) diff --git a/vignettes/von-bertalanffy.Rmd b/vignettes/von-bertalanffy.Rmd index 6b6824e..303bd54 100644 --- a/vignettes/von-bertalanffy.Rmd +++ b/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. From f53700ec44591ad5b3847dd6a7da5149fe6b3388 Mon Sep 17 00:00:00 2001 From: Tess-LaCoil Date: Wed, 4 Dec 2024 15:37:20 +1100 Subject: [PATCH 2/3] Updated namespace, documentation. --- NAMESPACE | 1 + man/hmde_plot_obs_est_inds.Rd | 44 +++++++++++++++++++++++++++++++++++ 2 files changed, 45 insertions(+) create mode 100644 man/hmde_plot_obs_est_inds.Rd diff --git a/NAMESPACE b/NAMESPACE index 73a94af..fa43685 100644 --- a/NAMESPACE +++ b/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) diff --git a/man/hmde_plot_obs_est_inds.Rd b/man/hmde_plot_obs_est_inds.Rd new file mode 100644 index 0000000..d359095 --- /dev/null +++ b/man/hmde_plot_obs_est_inds.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/hmde_plot_obs_est_inds.R +\name{hmde_plot_obs_est_inds} +\alias{hmde_plot_obs_est_inds} +\title{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.} +\usage{ +hmde_plot_obs_est_inds( + ind_id_vec = NULL, + n_ind_to_plot = NULL, + measurement_data = NULL, + xlab = "Time", + ylab = "Y(t)", + title = NULL +) +} +\arguments{ +\item{ind_id_vec}{vector with list of ind_id values} + +\item{n_ind_to_plot}{integer giving number of individuals to plot if not speciried} + +\item{measurement_data}{tibble with estimated measurements} + +\item{xlab}{character string for replacement x axis label} + +\item{ylab}{character string for replacement y axis label} + +\item{title}{character string for replacement plot title} + +\item{model}{model name character string} + +\item{colour}{character string for replacement line colour} + +\item{alpha}{real number for replacement alpha value} +} +\value{ +ggplot object +} +\description{ +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. +} From 04e3e00b7f433378b6f05fa37a92a064768f970a Mon Sep 17 00:00:00 2001 From: Tess-LaCoil Date: Wed, 4 Dec 2024 15:52:06 +1100 Subject: [PATCH 3/3] Removed unused args from hmde_plot_obs_est_ind function doc. --- R/hmde_plot_obs_est_inds.R | 3 --- man/hmde_plot_obs_est_inds.Rd | 6 ------ 2 files changed, 9 deletions(-) diff --git a/R/hmde_plot_obs_est_inds.R b/R/hmde_plot_obs_est_inds.R index 1e381a4..b1f0f65 100644 --- a/R/hmde_plot_obs_est_inds.R +++ b/R/hmde_plot_obs_est_inds.R @@ -2,15 +2,12 @@ #' on posterior estimates. Structured to take in the measurement_data tibble constructed by #' the hmde_extract_estimates function. #' -#' @param model model name character string #' @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 -#' @param colour character string for replacement line colour -#' @param alpha real number for replacement alpha value #' #' @return ggplot object #' @export diff --git a/man/hmde_plot_obs_est_inds.Rd b/man/hmde_plot_obs_est_inds.Rd index d359095..0940246 100644 --- a/man/hmde_plot_obs_est_inds.Rd +++ b/man/hmde_plot_obs_est_inds.Rd @@ -27,12 +27,6 @@ hmde_plot_obs_est_inds( \item{ylab}{character string for replacement y axis label} \item{title}{character string for replacement plot title} - -\item{model}{model name character string} - -\item{colour}{character string for replacement line colour} - -\item{alpha}{real number for replacement alpha value} } \value{ ggplot object