diff --git a/DESCRIPTION b/DESCRIPTION index 442f496..9a297bc 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -8,7 +8,7 @@ Authors@R: person(given = "Fonti", family = "Kar", role = c("ctb"), email = "f.kar@unsw.edu.au", comment = c(ORCID = "0000-0002-2760-3974")), person(given = "David", family= "Warton", role = c("aut", "ctb"), email = "david.warton@unsw.edu.au", comment = c(ORCID = "0000-0002-1642-628X")) ) -Description: What the package does (one paragraph). +Description: Wrapper for Stan that offers a number of in-built models to implement a hierarchical Bayesian longitudinal model for repeat observation data. Model choice selects the differential equation that is fit to the observations. Single and multi-individual models are available. License: MIT + file LICENSE Encoding: UTF-8 Roxygen: list(markdown = TRUE) @@ -18,13 +18,14 @@ Depends: R (>= 3.5.0) Imports: methods, + dplyr, + ggplot2, purrr, Rcpp (>= 0.12.0), RcppParallel (>= 5.0.1), rlang, rstan (>= 2.18.1), - rstantools (>= 2.3.1.1), - stats + rstantools (>= 2.3.1.1) LinkingTo: BH (>= 1.66.0), Rcpp (>= 0.12.0), @@ -35,7 +36,6 @@ LinkingTo: SystemRequirements: GNU make Suggests: knitr, - dplyr, rmarkdown, testthat (>= 3.0.0), withr, diff --git a/NAMESPACE b/NAMESPACE index 7ac7ab8..9167bb1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,11 +1,24 @@ # Generated by roxygen2: do not edit by hand 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_pars) +export(hmde_plot_de_pieces) +export(hmde_power_de) export(hmde_run) +export(hmde_vb_de) import(Rcpp) +import(dplyr) +import(ggplot2) import(methods) importFrom(RcppParallel,RcppParallelLibs) importFrom(rstan,sampling) importFrom(rstantools,rstan_config) +importFrom(stats,quantile) useDynLib(hmde, .registration = TRUE) diff --git a/R/hmde-package.R b/R/hmde-package.R index a39db48..7848f4d 100644 --- a/R/hmde-package.R +++ b/R/hmde-package.R @@ -8,6 +8,7 @@ #' @useDynLib hmde, .registration = TRUE #' @import methods #' @import Rcpp +#' @import dplyr #' @importFrom rstan sampling #' @importFrom rstantools rstan_config #' @importFrom RcppParallel RcppParallelLibs diff --git a/R/hmde_assign_data.R b/R/hmde_assign_data.R index 35c7150..c5d5a2e 100644 --- a/R/hmde_assign_data.R +++ b/R/hmde_assign_data.R @@ -1,30 +1,89 @@ -#' Assign data to template +#' Assign data to template for chosen model #' #' @param model_template output from hmde_model -#' @param ... data-masking name-value pairs +#' @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, ...){ - # Grab user expressions - user_code <- rlang::enquos(..., .check_assign = TRUE) +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.") + } + + if(!is.null(data)){ # Use provided tibble + user_fields <- names(data) + + } else { # Grab user expressions from individual list items and extract data + user_code <- rlang::enquos(..., .check_assign = TRUE) + user_fields <- names(user_code) + # Evaluate the RHS of expressions (the values) + data <- purrr::map(user_code, + ~rlang::eval_tidy(.x, env = rlang::caller_env()) + ) + } # Grab the names - fields <- names(user_code) + model_fields <- names(model_template) + + # Check user data has required names + if(grepl("multi", model_template$model)){ # Multi-individual with ind_id vec + for(i in c("ind_id", "time", "y_obs", "obs_index")){ + if(!i %in% user_fields){ + stop(paste0("Improper data structure: ", i, " missing")) + } + } - #TODO: Check names are in model_template + } else { # Single individual models + for(i in c("time", "y_obs", "obs_index")){ + if(!i %in% user_fields){ + stop(paste0("Improper data structure: ", i, " missing")) + } + } + } - # Evaluate the RHS of expressions (the values) - data <- purrr::map(user_code, - ~rlang::eval_tidy(.x, env = rlang::caller_env()) - ) + for(i in model_fields){ # Iterate through required fields and fill them + if(i %in% user_fields){ + model_template <- purrr::list_modify(model_template, !!!data[i]) + } else { + model_template[[i]] <- switch( + i, + step_size = step_size, + 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 + ) + } - for(i in fields){ - model_template <- purrr::list_modify(model_template, !!!data[i]) + if(is.null(model_template[[i]])){ #Report missing data + stop(paste("Improper data structure: Data missing:", i)) + } } - #TODO: Check if N is supplied, if not, give error + #Check lengths for y_obs, obs_index, time, ind_id + vec_lengths <- c( + model_template$n_obs, + length(model_template$y_obs), + length(model_template$obs_index), + length(model_template$time) + ) + if(grepl("multi", model_template$model)){ # Multi-individual with ind_id vector + vec_lengths[5] <- length(model_template$ind_id) + + #Check number ind ID values + ind_id_lengths <- c(model_template$n_ind,length(unique(data$ind_id))) + if(length(unique(ind_id_lengths))!=1){ + stop("Different values for n_ind and number of unique entries in ind_id.") + } + } + + if(length(unique(vec_lengths))!=1){ + stop("Improper data structure: Different lengths of data vectors.") + } return(model_template) } diff --git a/R/hmde_extract_estimates.R b/R/hmde_extract_estimates.R new file mode 100644 index 0000000..a088685 --- /dev/null +++ b/R/hmde_extract_estimates.R @@ -0,0 +1,175 @@ +#' 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 +#' +#' @return named list with data frames for measurement, individual, population-level, and error parameter estimates +#' @export +#' @import dplyr +#' @importFrom stats quantile + +hmde_extract_estimates <- function(model = NULL, + fit = NULL, + input_measurement_data = NULL){ + #Check for fit + if(is.null(fit)){ + stop("Fit not provided.") + } + + if(typeof(fit) != "S4"){ + stop("Fit not S4 stanfit type.") + } + + #Check for model + if(!model %in% hmde_model_name()){ + stop("Model name not recognised. Run hmde_model_name() to see available models.") + } + + #Check for input measurement data + for(i in c("y_obs", "time", "obs_index")){ + if(!i %in% names(input_measurement_data)){ + stop(paste("Input measurements information missing:", i)) + } + } + + + estimate_list <- list() + par_names <- hmde_model_pars(model) + + if(grepl("multi", model)){ #Get n_ind for multi-individual + n_ind <- length(unique(input_measurement_data$ind_id)) + } else { + n_ind <- 1 + } + + #Extract samples + samples <- rstan::extract(fit, permuted = TRUE, inc_warmup = FALSE) + sample_par_names <- names(samples) + + #Check parameter names + for(i in par_names[1:(length(par_names)-1)]){ + for(j in i){ + if(!j %in% sample_par_names){ + stop(paste( + "Parameter missing from model:", j + )) + } + } + } + + #Extract measurement, individual-level, and error parameter estimates and add to list + estimate_list$measurement_data <- hmde_extract_measurement_ests(samples, + par_names$measurement_pars_names, + input_measurement_data) + + estimate_list$individual_data <- hmde_extract_individual_par_ests(samples, + par_names$individual_pars_names, + n_ind) + + estimate_list$error_data <- hmde_extract_error_par_ests(samples, par_names$error_pars_names) + + #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) + } + + return(estimate_list) +} + + +#' Sample extraction for measurement-level estimates +#' @keywords internal +#' @noRd +hmde_extract_measurement_ests <- function(samples = NULL, + measurement_pars_names = NULL, + input_measurement_data = NULL){ + measurement_data <- input_measurement_data + + for(i in measurement_pars_names){ + measurement_data[[i]] <- apply(samples[[i]], 2, mean) + } + + return(measurement_data) +} + +#' Sample extraction for individual-level parameters +#' @keywords internal +#' @noRd +hmde_extract_individual_par_ests <- function(samples = NULL, + individual_pars_names = NULL, + n_ind = NULL){ + individual_data <- tibble(ind_id = 1:n_ind) + + #Extract mean of parameter posterior distributions + + if(n_ind > 1){ + for(i in individual_pars_names){ + individual_data[[paste0(i, "_mean")]] <- apply(samples[[i]], 2, mean) + individual_data[[paste0(i, "_median")]] <- apply(samples[[i]], 2, median) + individual_data[[paste0(i, "_CI_lower")]] <- apply(samples[[i]], 2, + stats::quantile, probs=c(0.025)) + individual_data[[paste0(i, "_CI_upper")]] <- apply(samples[[i]], 2, + stats::quantile, probs=c(0.975)) + } + } else { + for(i in individual_pars_names){ + individual_data[[paste0(i, "_mean")]] <- mean(samples[[i]]) + individual_data[[paste0(i, "_median")]] <- median(samples[[i]], 2, ) + individual_data[[paste0(i, "_CI_lower")]] <- as.numeric(stats::quantile(samples[[i]], + probs=c(0.025))) + individual_data[[paste0(i, "_CI_upper")]] <- as.numeric(stats::quantile(samples[[i]], + probs=c(0.975))) + } + } + + + return(individual_data) +} + +#' #' Sample extraction for population-level parameters +#' @keywords internal +#' @noRd +hmde_extract_pop_par_ests <- function(samples = NULL, + population_pars_names = NULL){ + population_data <- tibble() + + #Extract mean of parameter posterior distributions + for(i in population_pars_names){ + pop_data_temp <- tibble(par_name = i) + pop_data_temp$mean <- mean(samples[[i]]) + pop_data_temp$median <- median(samples[[i]]) + pop_data_temp$CI_lower <- as.numeric(stats::quantile(samples[[i]], + probs=c(0.025))) + pop_data_temp$CI_upper <- as.numeric(stats::quantile(samples[[i]], + probs=c(0.975))) + + population_data <- rbind(population_data, pop_data_temp) + } + + return(population_data) +} + +#' #' Sample extraction for error parameters +#' @keywords internal +#' @noRd +hmde_extract_error_par_ests <- function(samples = NULL, + error_pars_names = NULL){ + error_data <- tibble() + + #Extract mean of parameter posterior distributions + for(i in error_pars_names){ + error_data_temp <- tibble(par_name = i) + error_data_temp[["mean"]] <- mean(samples[[i]]) + error_data_temp[["median"]] <- median(samples[[i]]) + error_data_temp[["CI_lower"]] <- as.numeric(stats::quantile(samples[[i]], + probs=c(0.025))) + error_data_temp[["CI_upper"]] <- as.numeric(stats::quantile(samples[[i]], + probs=c(0.975))) + + error_data <- rbind(error_data, error_data_temp) + } + + return(error_data) +} diff --git a/R/hmde_model_des.R b/R/hmde_model_des.R new file mode 100644 index 0000000..7ee91d6 --- /dev/null +++ b/R/hmde_model_des.R @@ -0,0 +1,95 @@ +#' Function to select DE given model name +#' @param model character string model name +#' +#' @return DE function corresponding to specific model +#' @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.") + } + + output <- switch( + model, + constant_single_ind = hmde_const_de, + constant_multi_ind = hmde_const_de, + canham_single_ind = hmde_canham_de, + canham_multi_ind = hmde_canham_de, + power_single_ind = hmde_power_de, + power_multi_ind = hmde_power_de, + vb_single_ind = hmde_vb_de, + vb_multi_ind = hmde_vb_de, + linear_single_ind = hmde_linear_de + ) + + return(output) +} + +#' Differential equation for constant growth single and multi- individual models +#' @param y input real +#' @param pars list of parameter beta +#' +#' @return value of differential equation at y +#' @export + +hmde_const_de <- function(y = NULL, pars = NULL){ + return( + pars[[1]] + ) +} + + +#' Differential equation for Canham growth single and multi- individual models +#' @param y input real +#' @param pars list of parametera g_max, S_max, k +#' +#' @return value of differential equation at y +#' @export + +hmde_canham_de <- function(y = NULL, pars = NULL){ + return( + pars[[1]] * exp(-0.5 * + ((log(y / pars[[2]]) / pars[[3]]) ^ 2) + ) + ) +} + +#' Differential equation for power law growth single and multi- individual models +#' @param y input real +#' @param pars list of parameters coefficient, power +#' +#' @return value of differential equation at y +#' @export + +hmde_power_de <- function(y = NULL, pars = NULL){ + return( + pars[[1]] * (y ^ -pars[[2]]) + ) +} + + +#' Differential equation for von Bertalanffy growth single and multi- individual models +#' @param y input real +#' @param pars list of parameters beta, Y_max +#' +#' @return value of differential equation at y +#' @export + +hmde_vb_de <- function(y = NULL, pars = NULL){ + return( + pars[[2]] * (pars[[1]] - y) + ) +} + +#' Differential equation for linear 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){ + return( + pars[[1]] + pars[[2]] * y + ) +} diff --git a/R/hmde_model_names.R b/R/hmde_model_names.R new file mode 100644 index 0000000..c1e50b5 --- /dev/null +++ b/R/hmde_model_names.R @@ -0,0 +1,18 @@ +#' Returns names of available models. +#' +#' @return vector of character strings for model names. +#' @export + +hmde_model_name <- function(){ + output <- c("constant_single_ind", + "constant_multi_ind", + "canham_single_ind", + "canham_multi_ind", + "power_single_ind", + "power_multi_ind", + "vb_single_ind", + "vb_multi_ind", + "linear_single_ind") + + return(output) +} diff --git a/R/hmde_model_pars.R b/R/hmde_model_pars.R new file mode 100644 index 0000000..e1b3961 --- /dev/null +++ b/R/hmde_model_pars.R @@ -0,0 +1,133 @@ +#' Show parameter list for hmde supported model +#' +#' @param model model name character string +#' +#' @return named list that matches Stan model parameters +#' @export + +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.") + } + + output <- switch(model, + constant_single_ind = hmde_const_single_ind_pars(), + constant_multi_ind = hmde_const_multi_ind_pars(), + canham_single_ind = hmde_canham_single_ind_pars(), + canham_multi_ind = hmde_canham_multi_ind_pars(), + power_single_ind = hmde_power_single_ind_pars(), + power_multi_ind = hmde_power_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()) + + return(output) +} + +#' Parameter names for constant growth single individual model +#' @keywords internal +#' @noRd + +hmde_const_single_ind_pars <- function(){ + list(measurement_pars_names = c("y_hat", "Delta_hat"), + individual_pars_names = c("ind_beta"), + error_pars_names = c("global_error_sigma"), + model = "constant_single_ind") +} + +#' Parameter names for constant growth single pop model +#' @keywords internal +#' @noRd + +hmde_const_multi_ind_pars <- function(){ + list(measurement_pars_names = c("y_hat", "Delta_hat"), + individual_pars_names = c("ind_beta"), + population_pars_names = c("pop_beta_mu", "pop_beta_sigma"), + error_pars_names = c("global_error_sigma"), + model = "constant_multi_ind") +} + +#' Parameter names for Canham growth single individual model +#' @keywords internal +#' @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"), + error_pars_names = c("global_error_sigma"), + model = "canham_single_ind") +} + +#' Parameter names for Canham growth single pop model +#' @keywords internal +#' @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"), + population_pars_names = c("pop_max_growth_mu", "pop_max_growth_sigma", + "pop_diameter_at_max_growth_mu", "pop_diameter_at_max_growth_sigma", + "pop_k_mu", "pop_k_sigma"), + error_pars_names = c("global_error_sigma"), + model = "canham_multi_ind") +} + +#' Parameter names for power law growth single individual model +#' @keywords internal +#' @noRd + +hmde_power_single_ind_pars <- function(){ + list(measurement_pars_names = c("y_hat", "Delta_hat"), + individual_pars_names = c("ind_coeff", "ind_power"), + error_pars_names = c("global_error_sigma"), + model = "power_single_ind") +} + +#' Parameter names for power law growth single pop model +#' @keywords internal +#' @noRd + +hmde_power_multi_ind_pars <- function(){ + list(measurement_pars_names = c("y_hat", "Delta_hat"), + individual_pars_names = c("ind_coeff", "ind_power"), + population_pars_names = c("pop_coeff_mu", "pop_coeff_sigma", + "pop_power_mu", "pop_power_sigma"), + error_pars_names = c("global_error_sigma"), + model = "power_multi_ind") +} + +#' Parameter names for von Bertalanffy growth single individual model +#' @keywords internal +#' @noRd + +hmde_vb_single_ind_pars <- function(){ + list(measurement_pars_names = c("y_hat", "Delta_hat"), + individual_pars_names = c("ind_max_size", "ind_growth_rate"), + error_pars_names = c("global_error_sigma"), + model = "vb_single_ind") +} + +#' Parameter names for von Bertalanffy growth single pop model +#' @keywords internal +#' @noRd + +hmde_vb_multi_ind_pars <- function(){ + list(measurement_pars_names = c("y_hat", "Delta_hat"), + individual_pars_names = c("ind_max_size", "ind_growth_rate"), + population_pars_names = c("pop_max_size_mu", "pop_max_size_sigma", + "pop_growth_rate_mu", "pop_growth_rate_sigma"), + error_pars_names = c("global_error_sigma"), + model = "vb_multi_ind") +} + +#' Parameter names for linear growth single individual model +#' @keywords internal +#' @noRd +#' +hmde_linear_single_ind_pars <- function(){ + list(measurement_pars_names = c("y_hat", "Delta_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_models.R b/R/hmde_models.R index d642484..d6b7e63 100644 --- a/R/hmde_models.R +++ b/R/hmde_models.R @@ -7,7 +7,9 @@ hmde_model <- function(model=NULL){ - #TODO: Need a mechanism to check model requested in one that is supported by hmde + if(!model %in% hmde_model_name()){ + stop("Model name not recognised. Run hmde_model_name() to see available models.") + } output <- switch(model, constant_single_ind = hmde_const_single_ind(), @@ -25,17 +27,6 @@ hmde_model <- function(model=NULL){ return(output) } -#' Data configuration template for linear model -#' @keywords internal -#' @noRd - -hmde_lm <- function(){ - list(X = NULL, - Y = NULL, - N = NULL, - model = "linear") -} - #' Data configuration template for constant growth single individual model #' @keywords internal #' @noRd diff --git a/R/hmde_plot_de_pieces.R b/R/hmde_plot_de_pieces.R new file mode 100644 index 0000000..7e4970d --- /dev/null +++ b/R/hmde_plot_de_pieces.R @@ -0,0 +1,94 @@ +#' Plot pieces of chosen differential equation model for each individual. +#' Function piece will go from the first fitted size to the last. +#' Accepted ggplot arguments will change the axis labels, title, line colour, alpha +#' +#' @param model model name character string +#' @param individual_data tibble with estimated DE parameters +#' @param measurement_data tibble with estimated measurements +#' @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_de_pieces <- function(model = NULL, + individual_data = NULL, + measurement_data = NULL, + xlab = "Y(t)", + ylab = "f", + title = NULL, + colour = "#006600", + alpha = 0.2){ + #Check for model + if(!model %in% hmde_model_name()){ + stop("Model name not recognised. Run hmde_model_name() to see available models.") + } + + if(is.null(model)){ + stop("Model not provided.") + } + + if(is.null(individual_data)){ + stop("Individual parameter data not provided.") + } + + if(is.null(measurement_data)){ + stop("Measurement data not provided.") + } + + #Extract initial and final sizes for each individual + initial_and_final_vals <- measurement_data %>% + group_by(ind_id) %>% + arrange(obs_index) %>% + mutate(y_0 = first(y_hat), + y_final = last(y_hat)) %>% + ungroup() %>% + 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, + DE_function = hmde_model_des(model), + xlab = xlab, + ylab = ylab, + title = title, + colour = colour, + alpha = alpha) + + return(plot) +} + +#' Produce plot for plot_de_pieces +#' @keywords internal +#' @noRd +hmde_ggplot_de_pieces <- function(pars_data, + DE_function, + xlab, + ylab, + title, + colour, + alpha){ + plot <- ggplot() + + xlim(min(pars_data$y_0), max(pars_data$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")) + + for(i in 1:nrow(pars_data)){ + args_list <- list(pars=pars_data[i,]) + 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])) + } + + return(plot) +} diff --git a/R/hmde_run.R b/R/hmde_run.R index c5deb94..071d9d7 100644 --- a/R/hmde_run.R +++ b/R/hmde_run.R @@ -7,6 +7,10 @@ #' @export hmde_run <- function(model_template, ...) { + #Check for model + if(!model_template$model %in% hmde_model_name()){ + stop("Model name not recognised. Run hmde_model_name() to see available models.") + } # Detect model out <- switch(model_template$model, diff --git a/inst/stan/canham_multi_ind.stan b/inst/stan/canham_multi_ind.stan index 1f8bc9e..f00a2da 100644 --- a/inst/stan/canham_multi_ind.stan +++ b/inst/stan/canham_multi_ind.stan @@ -68,7 +68,7 @@ parameters { real ind_size_at_max_growth[n_ind]; real ind_k[n_ind]; - //Species level + //Population level real pop_max_growth_mean; real pop_max_growth_sd; real pop_size_at_max_growth_mean; diff --git a/inst/stan/constant_multi_ind.stan b/inst/stan/constant_multi_ind.stan index 126a091..16d0809 100644 --- a/inst/stan/constant_multi_ind.stan +++ b/inst/stan/constant_multi_ind.stan @@ -61,7 +61,7 @@ model { ind_beta ~ lognormal(pop_beta_mu, pop_beta_sigma); - //Species level + //Population level pop_beta_mu ~ normal(0.1, 1); pop_beta_sigma ~cauchy(0.1, 1); diff --git a/inst/stan/power_multi_ind.stan b/inst/stan/power_multi_ind.stan index 0d464ff..a6c9514 100644 --- a/inst/stan/power_multi_ind.stan +++ b/inst/stan/power_multi_ind.stan @@ -68,7 +68,7 @@ parameters { real ind_coeff[n_ind]; real ind_power[n_ind]; - //Species level + // Population level real pop_coeff_mean; real pop_coeff_sd; real pop_power_mean; @@ -110,7 +110,7 @@ model { ind_coeff ~lognormal(pop_coeff_mean, pop_coeff_sd); ind_power ~lognormal(pop_power_mean, pop_power_sd); - //Species level + // Population level pop_coeff_mean ~normal(0, 2); pop_coeff_sd ~cauchy(0, 2); pop_power_mean ~normal(0, 2); diff --git a/inst/stan/vb_multi_ind.stan b/inst/stan/vb_multi_ind.stan index b48b19e..f9da59e 100644 --- a/inst/stan/vb_multi_ind.stan +++ b/inst/stan/vb_multi_ind.stan @@ -67,7 +67,7 @@ parameters { real ind_growth_par[n_ind]; real ind_max_size[n_ind]; - //Species level + //Population level real pop_growth_par_mean; real pop_growth_par_sd; real pop_max_size_mean; @@ -109,7 +109,7 @@ model { ind_growth_par ~lognormal(pop_growth_par_mean, pop_growth_par_sd); ind_max_size ~lognormal(pop_max_size_mean, pop_max_size_sd); - //Species level + //Population level pop_growth_par_mean ~normal(0, 2); pop_growth_par_sd ~cauchy(0, 2); pop_max_size_mean ~normal(max(log(y_obs)), 2); diff --git a/man/hmde_assign_data.Rd b/man/hmde_assign_data.Rd index b05e853..c72e6de 100644 --- a/man/hmde_assign_data.Rd +++ b/man/hmde_assign_data.Rd @@ -2,18 +2,22 @@ % Please edit documentation in R/hmde_assign_data.R \name{hmde_assign_data} \alias{hmde_assign_data} -\title{Assign data to template} +\title{Assign data to template for chosen model} \usage{ -hmde_assign_data(model_template, ...) +hmde_assign_data(model_template, data = NULL, step_size = NULL, ...) } \arguments{ \item{model_template}{output from hmde_model} -\item{...}{data-masking name-value pairs} +\item{data}{Input data tibble with columns including time, y_obs, obs_index, and additionally ind_id for multi-individual models} + +\item{step_size}{Step size for numerical integration.} + +\item{...}{data-masking name-value pairs allowing specific input of elements} } \value{ updated named list with your data assigned to Stan model parameters } \description{ -Assign data to template +Assign data to template for chosen model } diff --git a/man/hmde_canham_de.Rd b/man/hmde_canham_de.Rd new file mode 100644 index 0000000..5c2d802 --- /dev/null +++ b/man/hmde_canham_de.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/hmde_model_des.R +\name{hmde_canham_de} +\alias{hmde_canham_de} +\title{Differential equation for Canham growth single and multi- individual models} +\usage{ +hmde_canham_de(y = NULL, pars = NULL) +} +\arguments{ +\item{y}{input real} + +\item{pars}{list of parametera g_max, S_max, k} +} +\value{ +value of differential equation at y +} +\description{ +Differential equation for Canham growth single and multi- individual models +} diff --git a/man/hmde_const_de.Rd b/man/hmde_const_de.Rd new file mode 100644 index 0000000..9709ef9 --- /dev/null +++ b/man/hmde_const_de.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/hmde_model_des.R +\name{hmde_const_de} +\alias{hmde_const_de} +\title{Differential equation for constant growth single and multi- individual models} +\usage{ +hmde_const_de(y = NULL, pars = NULL) +} +\arguments{ +\item{y}{input real} + +\item{pars}{list of parameter beta} +} +\value{ +value of differential equation at y +} +\description{ +Differential equation for constant growth single and multi- individual models +} diff --git a/man/hmde_extract_estimates.Rd b/man/hmde_extract_estimates.Rd new file mode 100644 index 0000000..93b89b4 --- /dev/null +++ b/man/hmde_extract_estimates.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/hmde_extract_estimates.R +\name{hmde_extract_estimates} +\alias{hmde_extract_estimates} +\title{Extract samples and return measurement, individual, and population-level estimates} +\usage{ +hmde_extract_estimates(model = NULL, fit = NULL, input_measurement_data = NULL) +} +\arguments{ +\item{model}{model name character string} + +\item{fit}{fitted model Stan fit} + +\item{input_measurement_data}{data used to fit the model with ind_id, y_obs, time, obs_index tibble} +} +\value{ +named list with data frames for measurement, individual, population-level, and error parameter estimates +} +\description{ +Extract samples and return measurement, individual, and population-level estimates +} diff --git a/man/hmde_linear_de.Rd b/man/hmde_linear_de.Rd new file mode 100644 index 0000000..560d84a --- /dev/null +++ b/man/hmde_linear_de.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/hmde_model_des.R +\name{hmde_linear_de} +\alias{hmde_linear_de} +\title{Differential equation for linear growth single individual model} +\usage{ +hmde_linear_de(y = NULL, pars = NULL) +} +\arguments{ +\item{y}{input real} + +\item{pars}{list of parameters beta_0, beta_1} +} +\value{ +value of differential equation at y +} +\description{ +Differential equation for linear growth single individual model +} diff --git a/man/hmde_model_des.Rd b/man/hmde_model_des.Rd new file mode 100644 index 0000000..cebc5e7 --- /dev/null +++ b/man/hmde_model_des.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/hmde_model_des.R +\name{hmde_model_des} +\alias{hmde_model_des} +\title{Function to select DE given model name} +\usage{ +hmde_model_des(model = NULL) +} +\arguments{ +\item{model}{character string model name} +} +\value{ +DE function corresponding to specific model +} +\description{ +Function to select DE given model name +} diff --git a/man/hmde_model_name.Rd b/man/hmde_model_name.Rd new file mode 100644 index 0000000..f8ecbb1 --- /dev/null +++ b/man/hmde_model_name.Rd @@ -0,0 +1,14 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/hmde_model_names.R +\name{hmde_model_name} +\alias{hmde_model_name} +\title{Returns names of available models.} +\usage{ +hmde_model_name() +} +\value{ +vector of character strings for model names. +} +\description{ +Returns names of available models. +} diff --git a/man/hmde_model_pars.Rd b/man/hmde_model_pars.Rd new file mode 100644 index 0000000..1282da9 --- /dev/null +++ b/man/hmde_model_pars.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/hmde_model_pars.R +\name{hmde_model_pars} +\alias{hmde_model_pars} +\title{Show parameter list for hmde supported model} +\usage{ +hmde_model_pars(model = NULL) +} +\arguments{ +\item{model}{model name character string} +} +\value{ +named list that matches Stan model parameters +} +\description{ +Show parameter list for hmde supported model +} diff --git a/man/hmde_plot_de_pieces.Rd b/man/hmde_plot_de_pieces.Rd new file mode 100644 index 0000000..ca5f2c8 --- /dev/null +++ b/man/hmde_plot_de_pieces.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/hmde_plot_de_pieces.R +\name{hmde_plot_de_pieces} +\alias{hmde_plot_de_pieces} +\title{Plot pieces of chosen differential equation model for each individual. +Function piece will go from the first fitted size to the last. +Accepted ggplot arguments will change the axis labels, title, line colour, alpha} +\usage{ +hmde_plot_de_pieces( + model = NULL, + individual_data = NULL, + measurement_data = NULL, + xlab = "Y(t)", + ylab = "f", + title = NULL, + colour = "#006600", + alpha = 0.2 +) +} +\arguments{ +\item{model}{model name character string} + +\item{individual_data}{tibble with estimated DE parameters} + +\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{colour}{character string for replacement line colour} + +\item{alpha}{real number for replacement alpha value} +} +\value{ +ggplot object +} +\description{ +Plot pieces of chosen differential equation model for each individual. +Function piece will go from the first fitted size to the last. +Accepted ggplot arguments will change the axis labels, title, line colour, alpha +} diff --git a/man/hmde_power_de.Rd b/man/hmde_power_de.Rd new file mode 100644 index 0000000..337f246 --- /dev/null +++ b/man/hmde_power_de.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/hmde_model_des.R +\name{hmde_power_de} +\alias{hmde_power_de} +\title{Differential equation for power law growth single and multi- individual models} +\usage{ +hmde_power_de(y = NULL, pars = NULL) +} +\arguments{ +\item{y}{input real} + +\item{pars}{list of parameters coefficient, power} +} +\value{ +value of differential equation at y +} +\description{ +Differential equation for power law growth single and multi- individual models +} diff --git a/man/hmde_vb_de.Rd b/man/hmde_vb_de.Rd new file mode 100644 index 0000000..160169b --- /dev/null +++ b/man/hmde_vb_de.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/hmde_model_des.R +\name{hmde_vb_de} +\alias{hmde_vb_de} +\title{Differential equation for von Bertalanffy growth single and multi- individual models} +\usage{ +hmde_vb_de(y = NULL, pars = NULL) +} +\arguments{ +\item{y}{input real} + +\item{pars}{list of parameters beta, Y_max} +} +\value{ +value of differential equation at y +} +\description{ +Differential equation for von Bertalanffy growth single and multi- individual models +} diff --git a/tests/testthat/test-hmde_assign_data.R b/tests/testthat/test-hmde_assign_data.R index 61761f2..93cd1b4 100644 --- a/tests/testthat/test-hmde_assign_data.R +++ b/tests/testthat/test-hmde_assign_data.R @@ -1,20 +1,123 @@ -test_that("Execution and output", { +test_that("Execution and output: Constant single ind", { constant_single <- hmde_model("constant_single_ind") |> hmde_assign_data(n_obs = 2, y_obs = c(1,2), obs_index = c(1,2), time = c(0,1), y_0_obs = 1) - + expect_named(constant_single) - + expect_visible(constant_single) - + expect_type(constant_single, "list") - + constant_empty <- hmde_model("constant_single_ind") - + expect_null(constant_empty$y_obs) - + expect_true(length(constant_single$y_obs) > 0) }) + +test_that("Execution and output: Constant multi ind manual input", { + #Full manual input + constant_multi <- hmde_model("constant_multi_ind") |> + hmde_assign_data(n_obs = nrow(Trout_Size_Data), #integer + n_ind = max(Trout_Size_Data$ind_id), #integer + y_obs = Trout_Size_Data$y_obs, #vector length N_obs + obs_index = Trout_Size_Data$obs_index, #vector length N_obs + time = Trout_Size_Data$time, #Vector length N_obs + ind_id = Trout_Size_Data$ind_id, #Vector length N_obs + y_0_obs = Trout_Size_Data$y_obs[which(Trout_Size_Data$obs_index == 1)] #Vector length N_ind + ) + + expect_named(constant_multi) + + expect_visible(constant_multi) + + expect_type(constant_multi, "list") + + constant_empty <- hmde_model("constant_multi_ind") + + expect_null(constant_empty$y_obs) + + expect_true(length(constant_multi$y_obs) > 0) +}) + +test_that("Execution and output: Constant multi ind minimum manual input", { + #Full manual input + constant_multi <- hmde_model("constant_multi_ind") |> + hmde_assign_data(y_obs = Trout_Size_Data$y_obs, #vector length N_obs + obs_index = Trout_Size_Data$obs_index, #vector length N_obs + time = Trout_Size_Data$time, #Vector length N_obs + ind_id = Trout_Size_Data$ind_id #Vector length N_obs + ) + + expect_named(constant_multi) + + expect_visible(constant_multi) + + expect_type(constant_multi, "list") + + expect_true(length(constant_multi$y_obs) > 0) +}) + +test_that("Execution and output: Constant multi ind tibble input", { + #Full manual input + constant_multi <- hmde_model("constant_multi_ind") |> + hmde_assign_data(data = Trout_Size_Data) + + expect_named(constant_multi) + + expect_visible(constant_multi) + + expect_type(constant_multi, "list") + + expect_true(length(constant_multi$y_obs) > 0) +}) + +test_that("Execution and output: bad input", { + #model does not exist + expect_error( + hmde_model("this_is_not_a_model") |> + hmde_assign_data(y_obs = Trout_Size_Data$y_obs, #vector length N_obs + obs_index = Trout_Size_Data$obs_index, #vector length N_obs + time = Trout_Size_Data$time, #Vector length N_obs + ind_id = Trout_Size_Data$ind_id, #Vector length N_obs + n_ind = 1 + ) + ) + + #Wrong number of individuals + expect_error( + hmde_model("constant_multi_ind") |> + hmde_assign_data(y_obs = Trout_Size_Data$y_obs, #vector length N_obs + obs_index = Trout_Size_Data$obs_index, #vector length N_obs + time = Trout_Size_Data$time, #Vector length N_obs + ind_id = Trout_Size_Data$ind_id, #Vector length N_obs + n_ind = 1 + ) + ) + + #Mismatched vector lengths + expect_error( + hmde_model("constant_multi_ind") |> + hmde_assign_data(y_obs = Trout_Size_Data$y_obs, #vector length N_obs + obs_index = Trout_Size_Data$obs_index[1:5], #vector length N_obs + time = Trout_Size_Data$time, #Vector length N_obs + ind_id = Trout_Size_Data$ind_id + ) + ) + + #Missing data + expect_error( + hmde_model("constant_multi_ind") |> + hmde_assign_data(data = c(0,1)) + ) + + #Wrong data set + expect_error( + hmde_model("constant_multi_ind") |> + hmde_assign_data(data = mtcars) + ) +}) diff --git a/tests/testthat/test-hmde_extract_estimates.R b/tests/testthat/test-hmde_extract_estimates.R new file mode 100644 index 0000000..f14202f --- /dev/null +++ b/tests/testthat/test-hmde_extract_estimates.R @@ -0,0 +1,25 @@ +test_that("Execution and output: extract_estimates function", { + data <- data.frame( + y_obs = c(1.1, 2.0, 2.9), + obs_index = c(1, 2, 3), + time = c(0, 1, 2), + ind_id = c(1, 1, 1) + ) + + suppressWarnings( + constant_single_fit <- hmde_model("constant_single_ind") |> + hmde_assign_data(data = data) |> + hmde_run(chains = 1, iter = 20, cores = 1, + verbose = FALSE, show_messages = FALSE) + ) + + output <- hmde_extract_estimates(model = "constant_single_ind", + fit = constant_single_fit, + input_measurement_data = data) + + expect_named(output) + + expect_visible(output) + + expect_type(output, "list") +}) diff --git a/tests/testthat/test-hmde_plot_de_pieces.R b/tests/testthat/test-hmde_plot_de_pieces.R new file mode 100644 index 0000000..9d7dfbd --- /dev/null +++ b/tests/testthat/test-hmde_plot_de_pieces.R @@ -0,0 +1,63 @@ +test_that("Execution and output: plot_de_pieces function", { + suppressWarnings( + multi_ind_trout <- hmde_model("constant_multi_ind") |> + hmde_assign_data(data = Trout_Size_Data) |> + hmde_run(chains = 1, iter = 50, cores = 1, + verbose = FALSE, show_messages = FALSE) + ) + + output <- hmde_extract_estimates(model = "constant_multi_ind", + fit = multi_ind_trout, + input_measurement_data = Trout_Size_Data) + + plot <- hmde_plot_de_pieces(model = "constant_multi_ind", + individual_data = output$individual_data, + measurement_data = output$measurement_data, + xlab = "S(t)", + ylab = "g", + title = "Constant growth") + + expect_named(plot) + + expect_visible(plot) + + expect_type(plot, "list") +}) + + +test_that("Execution and output: bad input", { + suppressWarnings( + multi_ind_trout <- hmde_model("constant_multi_ind") |> + hmde_assign_data(data = Trout_Size_Data) |> + hmde_run(chains = 1, iter = 50, cores = 1, + verbose = FALSE, show_messages = FALSE) + ) + + output <- hmde_extract_estimates(model = "constant_multi_ind", + fit = multi_ind_trout, + input_measurement_data = Trout_Size_Data) + + expect_error( + hmde_plot_de_pieces(model = "not_a_model", + individual_data = output$individual_data, + measurement_data = output$measurement_data) + ) + + expect_error( + hmde_plot_de_pieces(model = NULL, + individual_data = output$individual_data, + measurement_data = output$measurement_data) + ) + + expect_error( + hmde_plot_de_pieces(model = "constant_multi_ind", + individual_data = NULL, + measurement_data = output$measurement_data) + ) + + expect_error( + hmde_plot_de_pieces(model = "constant_multi_ind", + individual_data = output$individual_data, + measurement_data = NULL) + ) +})