diff --git a/NAMESPACE b/NAMESPACE index d5b46ca..4fe8aef 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,6 +4,7 @@ export(add_exposure_data) export(biokinetics) export(biokinetics_priors) export(convert_log2_scale_inverse) +export(plot_prior_predictive) importFrom(R6,R6Class) importFrom(data.table,":=") importFrom(data.table,.BY) diff --git a/R/biokinetics.R b/R/biokinetics.R index 6d4a558..b62a4e7 100644 --- a/R/biokinetics.R +++ b/R/biokinetics.R @@ -265,6 +265,17 @@ biokinetics <- R6::R6Class( package = "epikinetics" ) }, + #' @description Plot the kinetics trajectory predicted by the model priors. + #' @return A ggplot2 object. + #' @param tmax Integer. The number of time points in each simulated trajectory. Default 150. + #' @param n_draws Integer. The number of trajectories to simulate. Default 2000. + plot_prior_predictive = function(tmax = 150, + n_draws = 2000) { + plot_prior_predictive(private$priors, + tmax = tmax, + n_draws = n_draws, + data = private$data) + }, #' @description View the data that is passed to the stan model, for debugging purposes. #' @return A list of arguments that will be passed to the stan model. get_stan_data = function() { diff --git a/R/plot.R b/R/plot.R new file mode 100644 index 0000000..cc61250 --- /dev/null +++ b/R/plot.R @@ -0,0 +1,48 @@ +#' @title Simulate biomarker kinetics predicted by the given biokinetics priors +#' and optionally compare to a dataset. +#' @export +#' @description Simulate trajectories by drawing random samples from the given +#' priors for each parameter in the biokinetics model. +#' @return A ggplot2 object. +#' @param priors A named list of type 'biokinetics_priors'. +#' @param tmax Integer. The number of time points in each simulated trajectory. Default 150. +#' @param n_draws Integer. The number of trajectories to simulate. Default 2000. +#' @param data Optional data.frame with columns t_since_last_exp and value. The raw data to compare to. +plot_prior_predictive <- function(priors, + tmax = 150, + n_draws = 2000, + data = NULL) { + validate_priors(priors) + if (!is.null(data)) { + validate_required_cols(data, c("t_since_last_exp", "value")) + } + params <- data.table( + t0 = rnorm(n_draws, priors$mu_t0, priors$sigma_t0), # titre value at t0 + tp = rnorm(n_draws, priors$mu_tp, priors$sigma_tp), # time of peak + ts = rnorm(n_draws, priors$mu_ts, priors$sigma_ts), # time of set point + m1 = rnorm(n_draws, priors$mu_m1, priors$sigma_m1), # gradient 1 + m2 = rnorm(n_draws, priors$mu_m2, priors$sigma_m2), # gradient 2 + m3 = rnorm(n_draws, priors$mu_m3, priors$sigma_m3) # gradient 3 + ) + + times <- data.table(t = 1:tmax) + params_and_times <- times[, as.list(params), by = times] + + params_and_times[, mu := biokinetics_simulate_trajectory(t, t0, tp, ts, m1, m2, m3), + by = c("t", "t0", "tp", "ts", "m1", "m2", "m3")] + + summary <- params_and_times %>% + group_by(t) %>% + summarise(me = quantile(mu, 0.5, names = FALSE), + lo = quantile(mu, 0.025, names = FALSE), + hi = quantile(mu, 0.975, names = FALSE)) + + plot <- ggplot(summary) + + geom_line(aes(x = t, y = me)) + + geom_ribbon(aes(x = t, ymin = lo, ymax = hi), alpha = 0.5) + + if (!is.null(data)) { + plot <- plot + geom_point(data = data, aes(x = t_since_last_exp, y = value)) + } + plot +} diff --git a/man/biokinetics.Rd b/man/biokinetics.Rd index 5481497..9d6b3fd 100644 --- a/man/biokinetics.Rd +++ b/man/biokinetics.Rd @@ -11,6 +11,7 @@ fit it to a dataset. \subsection{Public methods}{ \itemize{ \item \href{#method-biokinetics-new}{\code{biokinetics$new()}} +\item \href{#method-biokinetics-plot_prior_predictive}{\code{biokinetics$plot_prior_predictive()}} \item \href{#method-biokinetics-get_stan_data}{\code{biokinetics$get_stan_data()}} \item \href{#method-biokinetics-get_covariate_lookup_table}{\code{biokinetics$get_covariate_lookup_table()}} \item \href{#method-biokinetics-fit}{\code{biokinetics$fit()}} @@ -62,6 +63,28 @@ An epikinetics::biokinetics object. } } \if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-biokinetics-plot_prior_predictive}{}}} +\subsection{Method \code{plot_prior_predictive()}}{ +Plot the kinetics trajectory predicted by the model priors. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{biokinetics$plot_prior_predictive(tmax = 150, n_draws = 2000)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{tmax}}{Integer. The number of time points in each simulated trajectory. Default 150.} + +\item{\code{n_draws}}{Integer. The number of trajectories to simulate. Default 2000.} +} +\if{html}{\out{
}} +} +\subsection{Returns}{ +A ggplot2 object. +} +} +\if{html}{\out{
}} \if{html}{\out{}} \if{latex}{\out{\hypertarget{method-biokinetics-get_stan_data}{}}} \subsection{Method \code{get_stan_data()}}{ diff --git a/man/plot_prior_predictive.Rd b/man/plot_prior_predictive.Rd new file mode 100644 index 0000000..c1f64f7 --- /dev/null +++ b/man/plot_prior_predictive.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot.R +\name{plot_prior_predictive} +\alias{plot_prior_predictive} +\title{Simulate biomarker kinetics predicted by the given biokinetics priors +and optionally compare to a dataset.} +\usage{ +plot_prior_predictive(priors, tmax = 150, n_draws = 2000, data = NULL) +} +\arguments{ +\item{priors}{A named list of type 'biokinetics_priors'.} + +\item{tmax}{Integer. The number of time points in each simulated trajectory. Default 150.} + +\item{n_draws}{Integer. The number of trajectories to simulate. Default 2000.} + +\item{data}{Optional data.frame with columns t_since_last_exp and value. The raw data to compare to.} +} +\value{ +A ggplot2 object. +} +\description{ +Simulate trajectories by drawing random samples from the given +priors for each parameter in the biokinetics model. +} diff --git a/tests/testthat.R b/tests/testthat.R index 9fc03dd..41aef3c 100644 --- a/tests/testthat.R +++ b/tests/testthat.R @@ -1,4 +1,5 @@ library(testthat) library(epikinetics) +library(vdiffr) test_check("epikinetics") diff --git a/tests/testthat/_snaps/plots/priorpredictive.svg b/tests/testthat/_snaps/plots/priorpredictive.svg new file mode 100644 index 0000000..44d0414 --- /dev/null +++ b/tests/testthat/_snaps/plots/priorpredictive.svg @@ -0,0 +1,2310 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +0 +5 +10 + + + + + + + +0 +200 +400 +600 +t +me +priorpredictive + + diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R new file mode 100644 index 0000000..3db51f8 --- /dev/null +++ b/tests/testthat/test-plots.R @@ -0,0 +1,41 @@ +test_that("Can plot prior prediction up to tmax", { + priors <- biokinetics_priors() + plot <- plot_prior_predictive(priors, tmax = 100, n_draws = 500) + expect_equal(nrow(plot$data), 100) + expect_equal(length(plot$layers), 2) +}) + +test_that("Can plot prior prediction with data points", { + data <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) + priors <- biokinetics_priors() + expect_error(plot_prior_predictive(priors, data = data), "Missing required columns: t_since_last_exp") + data[, `:=`(t_since_last_exp = as.integer(day - last_exp_day, units = "days"))] + plot <- plot_prior_predictive(priors, data = data, n_draws = 500) + expect_equal(length(plot$layers), 3) +}) + +test_that("Can plot prior predictions from model", { + data <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) + priors <- biokinetics_priors(mu_values = c(4.1, 11, 65, 0.2, -0.01, 0.01), + sigma_values = c(2.0, 2.0, 3.0, 0.01, 0.01, 0.001)) + + mod <- biokinetics$new(priors = priors, + data = data) + set.seed(1) + plot <- mod$plot_prior_predictive(tmax = 400, n_draws = 500) + expect_equal(nrow(plot$data), 400) + expect_equal(length(plot$layers), 3) +}) + +test_that("Prior predictions from model are the same", { + skip_on_ci() + data <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) + priors <- biokinetics_priors(mu_values = c(4.1, 11, 65, 0.2, -0.01, 0.01), + sigma_values = c(2.0, 2.0, 3.0, 0.01, 0.01, 0.001)) + + mod <- biokinetics$new(priors = priors, + data = data) + set.seed(1) + plot <- mod$plot_prior_predictive(tmax = 400, n_draws = 500) + vdiffr::expect_doppelganger("priorpredictive", plot) +})