From 730d6c4ada63cde9019c3ae5f7a8a3e690e4ab3e Mon Sep 17 00:00:00 2001 From: Tess O'Brien Date: Thu, 5 Dec 2024 16:45:02 +1100 Subject: [PATCH] =?UTF-8?q?Added=20a=20function=20to=20plot=20sizes=20over?= =?UTF-8?q?=20time=20in=20parallel=20with=20the=20DE=20plot=E2=80=A6=20(#4?= =?UTF-8?q?8)=20(#49)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Added a function to plot sizes over time in parallel with the DE plot function. Changed wording in README to be more accessible. * Updated namespace, documentation. * Removed unused args from hmde_plot_obs_est_ind function doc. --- NAMESPACE | 1 + R/hmde_plot_obs_est_inds.R | 78 ++++++++++++++++++++ README.md | 3 +- man/hmde_plot_obs_est_inds.Rd | 38 ++++++++++ tests/testthat/test-hmde_plot_obs_est_inds.R | 30 ++++++++ vignettes/canham.Rmd | 17 +---- vignettes/constant-growth.Rmd | 18 +---- vignettes/von-bertalanffy.Rmd | 17 +---- 8 files changed, 155 insertions(+), 47 deletions(-) create mode 100644 R/hmde_plot_obs_est_inds.R create mode 100644 man/hmde_plot_obs_est_inds.Rd create mode 100644 tests/testthat/test-hmde_plot_obs_est_inds.R diff --git a/NAMESPACE b/NAMESPACE index 2898f1e..733f919 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,6 +10,7 @@ export(hmde_model_des) 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) diff --git a/R/hmde_plot_obs_est_inds.R b/R/hmde_plot_obs_est_inds.R new file mode 100644 index 0000000..b1f0f65 --- /dev/null +++ b/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) +} 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/man/hmde_plot_obs_est_inds.Rd b/man/hmde_plot_obs_est_inds.Rd new file mode 100644 index 0000000..0940246 --- /dev/null +++ b/man/hmde_plot_obs_est_inds.Rd @@ -0,0 +1,38 @@ +% 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} +} +\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. +} 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 3168e3f..9ad321f 100644 --- a/vignettes/constant-growth.Rmd +++ b/vignettes/constant-growth.Rmd @@ -299,22 +299,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 bf1d2d4..5815f41 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.