From a4e479a0cae4c09d630fd8c64aaf08a21d511717 Mon Sep 17 00:00:00 2001 From: "alex.hill@gmail.com" Date: Mon, 21 Oct 2024 19:42:24 +0100 Subject: [PATCH 01/12] add prior predictive plot --- NAMESPACE | 1 + R/biokinetics.R | 11 + R/plot.R | 48 + man/biokinetics.Rd | 23 + man/plot_prior_predictive.Rd | 25 + tests/testthat.R | 1 + .../testthat/_snaps/plots/priorpredictive.svg | 2310 +++++++++++++++++ tests/testthat/test-plots.R | 41 + 8 files changed, 2460 insertions(+) create mode 100644 R/plot.R create mode 100644 man/plot_prior_predictive.Rd create mode 100644 tests/testthat/_snaps/plots/priorpredictive.svg create mode 100644 tests/testthat/test-plots.R 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) +}) From da15415cf1fc0978b4be2099fad92fb2dbcbf618 Mon Sep 17 00:00:00 2001 From: "alex.hill@gmail.com" Date: Mon, 21 Oct 2024 20:39:16 +0100 Subject: [PATCH 02/12] add input data plot imports from ggplot2 --- DESCRIPTION | 3 +- NAMESPACE | 9 + R/biokinetics.R | 5 + R/epikinetics-package.R | 1 + R/plot.R | 17 +- man/biokinetics.Rd | 14 + tests/testthat/_snaps/plots/inputdata.svg | 2408 +++++++++++++++++++++ tests/testthat/test-plots.R | 8 + 8 files changed, 2459 insertions(+), 6 deletions(-) create mode 100644 tests/testthat/_snaps/plots/inputdata.svg diff --git a/DESCRIPTION b/DESCRIPTION index 90156c1..f7d12c9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -13,6 +13,7 @@ Imports: data.table, forcats, fs, + ggplot2, instantiate, logger, mosaic, @@ -26,11 +27,11 @@ Additional_repositories: SystemRequirements: CmdStan (https://mc-stan.org/users/interfaces/cmdstan) Suggests: dplyr, - ggplot2, knitr, lubridate, rmarkdown, testthat + vdiffr VignetteBuilder: knitr LinkingTo: cpp11 diff --git a/NAMESPACE b/NAMESPACE index 4fe8aef..4a255a5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -15,4 +15,13 @@ importFrom(data.table,.N) importFrom(data.table,.NGRP) importFrom(data.table,.SD) importFrom(data.table,data.table) +importFrom(ggplot2,aes) +importFrom(ggplot2,facet_wrap) +importFrom(ggplot2,geom_line) +importFrom(ggplot2,geom_point) +importFrom(ggplot2,geom_ribbon) +importFrom(ggplot2,geom_smooth) +importFrom(ggplot2,ggplot) +importFrom(ggplot2,guide_legend) +importFrom(ggplot2,guides) useDynLib(epikinetics, .registration = TRUE) diff --git a/R/biokinetics.R b/R/biokinetics.R index b62a4e7..eb3f7f3 100644 --- a/R/biokinetics.R +++ b/R/biokinetics.R @@ -276,6 +276,11 @@ biokinetics <- R6::R6Class( n_draws = n_draws, data = private$data) }, + #' @description Plot model input data with a smoothing function. + #' @return A ggplot2 object. + plot_data = function() { + plot_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/epikinetics-package.R b/R/epikinetics-package.R index 5ac6050..4404fc3 100644 --- a/R/epikinetics-package.R +++ b/R/epikinetics-package.R @@ -11,6 +11,7 @@ #' @importFrom data.table .NGRP #' @importFrom data.table .SD #' @importFrom data.table data.table +#' @importFrom ggplot2 aes facet_wrap geom_point geom_ribbon geom_line geom_smooth ggplot guides guide_legend #' @useDynLib epikinetics, .registration = TRUE ## usethis namespace: end diff --git a/R/plot.R b/R/plot.R index cc61250..3ebc65f 100644 --- a/R/plot.R +++ b/R/plot.R @@ -31,11 +31,9 @@ plot_prior_predictive <- function(priors, 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)) + summary <- params_and_times[, .(me = stats::quantile(mu, 0.5, names = FALSE), + lo = stats::quantile(mu, 0.025, names = FALSE), + hi = stats::quantile(mu, 0.975, names = FALSE)), by = t] plot <- ggplot(summary) + geom_line(aes(x = t, y = me)) + @@ -46,3 +44,12 @@ plot_prior_predictive <- function(priors, } plot } + +plot_data <- function(data) { + validate_required_cols(data, c("t_since_last_exp", "value", "titre_type")) + ggplot(data) + + geom_point(aes(x = t_since_last_exp, y = value, colour = titre_type)) + + geom_smooth(aes(x = t_since_last_exp, y = value, colour = titre_type)) + + facet_wrap(~titre_type) + + guides(colour = guide_legend(title = "Titre type")) +} diff --git a/man/biokinetics.Rd b/man/biokinetics.Rd index 9d6b3fd..842b49f 100644 --- a/man/biokinetics.Rd +++ b/man/biokinetics.Rd @@ -12,6 +12,7 @@ fit it to a dataset. \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-plot_data}{\code{biokinetics$plot_data()}} \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()}} @@ -80,6 +81,19 @@ Plot the kinetics trajectory predicted by the model priors. } \if{html}{\out{}} } +\subsection{Returns}{ +A ggplot2 object. +} +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-biokinetics-plot_data}{}}} +\subsection{Method \code{plot_data()}}{ +Plot model input data with a smoothing function. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{biokinetics$plot_data()}\if{html}{\out{
}} +} + \subsection{Returns}{ A ggplot2 object. } diff --git a/tests/testthat/_snaps/plots/inputdata.svg b/tests/testthat/_snaps/plots/inputdata.svg new file mode 100644 index 0000000..c0b05fd --- /dev/null +++ b/tests/testthat/_snaps/plots/inputdata.svg @@ -0,0 +1,2408 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Alpha + + + + + + + + + + +Ancestral + + + + + + + + + + +Delta + + + + + + +0 +200 +400 +600 + + + + +0 +200 +400 +600 + + + + +0 +200 +400 +600 +0 +3 +6 +9 + + + + +t_since_last_exp +value + +Titre type + + + + + + + + + + + + +Alpha +Ancestral +Delta +inputdata + + diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R index 3db51f8..d01a8e1 100644 --- a/tests/testthat/test-plots.R +++ b/tests/testthat/test-plots.R @@ -39,3 +39,11 @@ test_that("Prior predictions from model are the same", { plot <- mod$plot_prior_predictive(tmax = 400, n_draws = 500) vdiffr::expect_doppelganger("priorpredictive", plot) }) + +test_that("Plotted data is are the same", { + skip_on_ci() + data <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) + mod <- biokinetics$new(data = data) + plot <- mod$plot_data() + vdiffr::expect_doppelganger("inputdata", plot) +}) From d137fb217fdb7b582e6101c827ab00d8b5feb5f7 Mon Sep 17 00:00:00 2001 From: "alex.hill@gmail.com" Date: Tue, 22 Oct 2024 13:54:17 +0100 Subject: [PATCH 03/12] trajectory plots and doppelganger tests --- DESCRIPTION | 2 +- NAMESPACE | 2 + R/biokinetics.R | 35 +- R/epikinetics-package.R | 2 +- R/plot.R | 84 +- R/utils.R | 2 +- man/biokinetics.Rd | 16 +- ...lot.biokinetics_population_trajectories.Rd | 19 + man/plot_prior_predictive.Rd | 2 +- .../_snaps/plots/inputdata-covariates.svg | 2506 +++++++++++++++++ tests/testthat/_snaps/plots/inputdata.svg | 2 +- .../plots/populationtrajectories-data.svg | 2500 ++++++++++++++++ .../_snaps/plots/populationtrajectories.svg | 245 ++ tests/testthat/test-data.R | 4 +- tests/testthat/test-plots.R | 46 +- 15 files changed, 5420 insertions(+), 47 deletions(-) create mode 100644 man/plot.biokinetics_population_trajectories.Rd create mode 100644 tests/testthat/_snaps/plots/inputdata-covariates.svg create mode 100644 tests/testthat/_snaps/plots/populationtrajectories-data.svg create mode 100644 tests/testthat/_snaps/plots/populationtrajectories.svg diff --git a/DESCRIPTION b/DESCRIPTION index f7d12c9..bcab1f4 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -30,7 +30,7 @@ Suggests: knitr, lubridate, rmarkdown, - testthat + testthat, vdiffr VignetteBuilder: knitr LinkingTo: diff --git a/NAMESPACE b/NAMESPACE index 4a255a5..7fa16cf 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,6 @@ # Generated by roxygen2: do not edit by hand +S3method(plot,biokinetics_population_trajectories) export(add_exposure_data) export(biokinetics) export(biokinetics_priors) @@ -24,4 +25,5 @@ importFrom(ggplot2,geom_smooth) importFrom(ggplot2,ggplot) importFrom(ggplot2,guide_legend) importFrom(ggplot2,guides) +importFrom(ggplot2,scale_y_continuous) useDynLib(epikinetics, .registration = TRUE) diff --git a/R/biokinetics.R b/R/biokinetics.R index eb3f7f3..3f6d48a 100644 --- a/R/biokinetics.R +++ b/R/biokinetics.R @@ -142,6 +142,11 @@ biokinetics <- R6::R6Class( dt_out[, t_id := NULL] + if (private$scale == "natural") { + dt_out <- convert_log2_scale_inverse( + dt_out, vars_to_transform = "mu") + } + if (summarise == TRUE) { logger::log_info("Summarising into quantiles") dt_out <- summarise_draws( @@ -266,6 +271,8 @@ biokinetics <- R6::R6Class( ) }, #' @description Plot the kinetics trajectory predicted by the model priors. + #' Note that this is on a log scale, regardless of whether the data was provided + #' on a log or a natural scale. #' @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. @@ -276,10 +283,12 @@ biokinetics <- R6::R6Class( n_draws = n_draws, data = private$data) }, - #' @description Plot model input data with a smoothing function. + #' @description Plot model input data with a smoothing function. Note that + #' this plot is on a log scale, regardless of whether data was provided on a + #' log or a natural scale. #' @return A ggplot2 object. - plot_data = function() { - plot_data(private$data) + plot_model_inputs = function() { + plot_data(private$data, private$all_formula_vars) }, #' @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. @@ -388,6 +397,7 @@ biokinetics <- R6::R6Class( t_max = 150, summarise = TRUE, n_draws = 2500) { + private$check_fitted() validate_numeric(t_max) validate_logical(summarise) @@ -404,17 +414,9 @@ biokinetics <- R6::R6Class( dt_out <- dt_out[ , lapply(.SD, function(x) if (is.factor(x)) forcats::fct_drop(x) else x)] - if (private$scale == "log") { - return(dt_out) - } - - if (summarise) { - dt_out <- convert_log2_scale_inverse( - dt_out, vars_to_transform = c("me", "lo", "hi")) - } else { - dt_out <- convert_log2_scale_inverse( - dt_out, vars_to_transform = "mu") - } + class(dt_out) <- append("biokinetics_population_trajectories", class(dt_out)) + attr(dt_out, "summarised") <- summarise + attr(dt_out, "covariates") <- private$all_formula_vars dt_out }, #' @description Process the stan model results into a data.table. @@ -550,7 +552,10 @@ biokinetics <- R6::R6Class( by = c("calendar_day", "titre_type")) } - dt_out[, time_shift := time_shift] + dt_out <- dt_out[, time_shift := time_shift] + class(dt_out) <- append("biokinetics_individual_trajectories", class(dt_out)) + attr(dt_out, "summarised") <- summarise + dt_out } ) ) diff --git a/R/epikinetics-package.R b/R/epikinetics-package.R index 4404fc3..35b5299 100644 --- a/R/epikinetics-package.R +++ b/R/epikinetics-package.R @@ -11,7 +11,7 @@ #' @importFrom data.table .NGRP #' @importFrom data.table .SD #' @importFrom data.table data.table -#' @importFrom ggplot2 aes facet_wrap geom_point geom_ribbon geom_line geom_smooth ggplot guides guide_legend +#' @importFrom ggplot2 aes facet_wrap geom_point geom_ribbon geom_line geom_smooth ggplot guides guide_legend scale_y_continuous #' @useDynLib epikinetics, .registration = TRUE ## usethis namespace: end diff --git a/R/plot.R b/R/plot.R index 3ebc65f..ed11e5b 100644 --- a/R/plot.R +++ b/R/plot.R @@ -7,22 +7,28 @@ #' @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. +#' @param data Optional data.frame with columns time_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) + + # Declare variables to suppress notes when compiling package + # https://github.com/Rdatatable/data.table/issues/850#issuecomment-259466153 + t0 <- tp <- ts <- m1 <- m2 <- m3 <- NULL + time_since_last_exp <- me <- lo <- hi <- value <- mu <- NULL + if (!is.null(data)) { - validate_required_cols(data, c("t_since_last_exp", "value")) + validate_required_cols(data, c("time_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 + t0 = stats::rnorm(n_draws, priors$mu_t0, priors$sigma_t0), # titre value at t0 + tp = stats::rnorm(n_draws, priors$mu_tp, priors$sigma_tp), # time of peak + ts = stats::rnorm(n_draws, priors$mu_ts, priors$sigma_ts), # time of set point + m1 = stats::rnorm(n_draws, priors$mu_m1, priors$sigma_m1), # gradient 1 + m2 = stats::rnorm(n_draws, priors$mu_m2, priors$sigma_m2), # gradient 2 + m3 = stats::rnorm(n_draws, priors$mu_m3, priors$sigma_m3) # gradient 3 ) times <- data.table(t = 1:tmax) @@ -40,16 +46,66 @@ plot_prior_predictive <- function(priors, 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 <- plot + geom_point(data = data, aes(x = time_since_last_exp, y = value)) } plot } -plot_data <- function(data) { - validate_required_cols(data, c("t_since_last_exp", "value", "titre_type")) +plot_data <- function(data, covariates) { + validate_required_cols(data, c("time_since_last_exp", "value", "titre_type")) + # Declare variables to suppress notes when compiling package + # https://github.com/Rdatatable/data.table/issues/850#issuecomment-259466153 + time_since_last_exp <- value <- titre_type <- NULL + ggplot(data) + - geom_point(aes(x = t_since_last_exp, y = value, colour = titre_type)) + - geom_smooth(aes(x = t_since_last_exp, y = value, colour = titre_type)) + - facet_wrap(~titre_type) + + geom_point(aes(x = time_since_last_exp, y = value, colour = titre_type)) + + geom_smooth(aes(x = time_since_last_exp, y = value, colour = titre_type)) + + facet_wrap(eval(parse(text = facet_formula(covariates)))) + guides(colour = guide_legend(title = "Titre type")) } + +#' Plot method for "biokinetics_population_trajectories" class +#' +#' @param x An object of class "biokinetics_population_trajectories". These are +#' generated by running biokinetics$simulate_populate_trajectories(). See +#' \href{../../epikinetics/html/biokinetics.html#method-biokinetics-simulate_population_trajectories}{\code{biokinetics$simulate_population_trajectories()}} +#' @param \dots Further arguments passed to the method. For now the only extra argument supported is "data", which should be +#' a data.table containing raw data as provided to the biokinetics model. +#' @export +plot.biokinetics_population_trajectories <- function(x, ...) { + if (!attr(x, "summarised")) { + by <- setdiff(colnames(x), c("t0_pop", "tp_pop", "ts_pop", + "m1_pop", "m2_pop", "m3_pop", + ".draw", "mu")) + x <- summarise_draws( + x, column_name = "mu", by = by) + } + args <- list(...) + # Declare variables to suppress notes when compiling package + # https://github.com/Rdatatable/data.table/issues/850#issuecomment-259466153 + time_since_last_exp <- value <- me <- titre_type <- lo <- hi <- NULL + day <- last_exp_day <- NULL + + plot <- ggplot(x) + + geom_line(aes(x = time_since_last_exp, y = me, colour = titre_type)) + + geom_ribbon(aes(x = time_since_last_exp, + ymin = lo, + ymax = hi, + fill = titre_type), alpha = 0.5) + + if (!is.null(args$data)) { + data <- args$data + validate_required_cols(data) + plot <- plot + + geom_point(data = data, + aes(x = as.integer(day - last_exp_day, units = "days"), + y = value), size = 0.5, alpha = 0.5) + } + plot + scale_y_continuous(trans = "log2") + + facet_wrap(eval(parse(text = facet_formula(attr(x, "covariates"))))) + + guides(fill = guide_legend(title = "Titre type"), colour = "none") +} + +facet_formula <- function(covariates) { + paste("~", paste(c("titre_type", covariates), collapse = "+")) +} diff --git a/R/utils.R b/R/utils.R index 1ee1a2e..6d770f4 100644 --- a/R/utils.R +++ b/R/utils.R @@ -1,5 +1,5 @@ convert_log2_scale <- function( - dt_in, vars_to_transform = "titre", + dt_in, vars_to_transform = "value", simplify_limits = TRUE) { dt_out <- data.table::copy(dt_in) diff --git a/man/biokinetics.Rd b/man/biokinetics.Rd index 842b49f..feb9501 100644 --- a/man/biokinetics.Rd +++ b/man/biokinetics.Rd @@ -12,7 +12,7 @@ fit it to a dataset. \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-plot_data}{\code{biokinetics$plot_data()}} +\item \href{#method-biokinetics-plot_model_inputs}{\code{biokinetics$plot_model_inputs()}} \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()}} @@ -68,6 +68,8 @@ An epikinetics::biokinetics object. \if{latex}{\out{\hypertarget{method-biokinetics-plot_prior_predictive}{}}} \subsection{Method \code{plot_prior_predictive()}}{ Plot the kinetics trajectory predicted by the model priors. +Note that this is on a log scale, regardless of whether the data was provided +on a log or a natural scale. \subsection{Usage}{ \if{html}{\out{
}}\preformatted{biokinetics$plot_prior_predictive(tmax = 150, n_draws = 2000)}\if{html}{\out{
}} } @@ -86,12 +88,14 @@ A ggplot2 object. } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-biokinetics-plot_data}{}}} -\subsection{Method \code{plot_data()}}{ -Plot model input data with a smoothing function. +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-biokinetics-plot_model_inputs}{}}} +\subsection{Method \code{plot_model_inputs()}}{ +Plot model input data with a smoothing function. Note that +this plot is on a log scale, regardless of whether data was provided on a +log or a natural scale. \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{biokinetics$plot_data()}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{biokinetics$plot_model_inputs()}\if{html}{\out{
}} } \subsection{Returns}{ diff --git a/man/plot.biokinetics_population_trajectories.Rd b/man/plot.biokinetics_population_trajectories.Rd new file mode 100644 index 0000000..aeecc78 --- /dev/null +++ b/man/plot.biokinetics_population_trajectories.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot.R +\name{plot.biokinetics_population_trajectories} +\alias{plot.biokinetics_population_trajectories} +\title{Plot method for "biokinetics_population_trajectories" class} +\usage{ +\method{plot}{biokinetics_population_trajectories}(x, ...) +} +\arguments{ +\item{x}{An object of class "biokinetics_population_trajectories". These are +generated by running biokinetics$simulate_populate_trajectories(). See +\href{../../epikinetics/html/biokinetics.html#method-biokinetics-simulate_population_trajectories}{\code{biokinetics$simulate_population_trajectories()}}} + +\item{\dots}{Further arguments passed to the method. For now the only extra argument supported is "data", which should be +a data.table containing raw data as provided to the biokinetics model.} +} +\description{ +Plot method for "biokinetics_population_trajectories" class +} diff --git a/man/plot_prior_predictive.Rd b/man/plot_prior_predictive.Rd index c1f64f7..fc6aa0c 100644 --- a/man/plot_prior_predictive.Rd +++ b/man/plot_prior_predictive.Rd @@ -14,7 +14,7 @@ plot_prior_predictive(priors, tmax = 150, n_draws = 2000, data = NULL) \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.} +\item{data}{Optional data.frame with columns time_since_last_exp and value. The raw data to compare to.} } \value{ A ggplot2 object. diff --git a/tests/testthat/_snaps/plots/inputdata-covariates.svg b/tests/testthat/_snaps/plots/inputdata-covariates.svg new file mode 100644 index 0000000..20b5d78 --- /dev/null +++ b/tests/testthat/_snaps/plots/inputdata-covariates.svg @@ -0,0 +1,2506 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Ancestral + +Previously infected (Pre-Omicron) + + + + + + + + + + +Delta + +Infection naive + + + + + + + + + + +Delta + +Previously infected (Pre-Omicron) + + + + + + + + + + +Alpha + +Infection naive + + + + + + + + + + +Alpha + +Previously infected (Pre-Omicron) + + + + + + + + + + +Ancestral + +Infection naive + + + + + + +0 +200 +400 +600 + + + + +0 +200 +400 +600 + + + + +0 +200 +400 +600 +0 +3 +6 +9 + + + + +0 +3 +6 +9 + + + + +time_since_last_exp +value + +Titre type + + + + + + + + + + + + +Alpha +Ancestral +Delta +inputdata_covariates + + diff --git a/tests/testthat/_snaps/plots/inputdata.svg b/tests/testthat/_snaps/plots/inputdata.svg index c0b05fd..e4c509f 100644 --- a/tests/testthat/_snaps/plots/inputdata.svg +++ b/tests/testthat/_snaps/plots/inputdata.svg @@ -2384,7 +2384,7 @@ -t_since_last_exp +time_since_last_exp value Titre type diff --git a/tests/testthat/_snaps/plots/populationtrajectories-data.svg b/tests/testthat/_snaps/plots/populationtrajectories-data.svg new file mode 100644 index 0000000..6340bda --- /dev/null +++ b/tests/testthat/_snaps/plots/populationtrajectories-data.svg @@ -0,0 +1,2500 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Ancestral + +Previously infected (Pre-Omicron) + + + + + + + + + + +Delta + +Infection naive + + + + + + + + + + +Delta + +Previously infected (Pre-Omicron) + + + + + + + + + + +Alpha + +Infection naive + + + + + + + + + + +Alpha + +Previously infected (Pre-Omicron) + + + + + + + + + + +Ancestral + +Infection naive + + + + + + +0 +200 +400 +600 + + + + +0 +200 +400 +600 + + + + +0 +200 +400 +600 +128 +8192 +524288 +33554432 + + + + +128 +8192 +524288 +33554432 + + + + +time_since_last_exp +me + +Titre type + + + + + + +Alpha +Ancestral +Delta +populationtrajectories_data + + diff --git a/tests/testthat/_snaps/plots/populationtrajectories.svg b/tests/testthat/_snaps/plots/populationtrajectories.svg new file mode 100644 index 0000000..3c2d3b9 --- /dev/null +++ b/tests/testthat/_snaps/plots/populationtrajectories.svg @@ -0,0 +1,245 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Ancestral + +Previously infected (Pre-Omicron) + + + + + + + + + + +Delta + +Infection naive + + + + + + + + + + +Delta + +Previously infected (Pre-Omicron) + + + + + + + + + + +Alpha + +Infection naive + + + + + + + + + + +Alpha + +Previously infected (Pre-Omicron) + + + + + + + + + + +Ancestral + +Infection naive + + + + + + +0 +50 +100 +150 + + + + +0 +50 +100 +150 + + + + +0 +50 +100 +150 +128 +8192 +524288 +33554432 + + + + +128 +8192 +524288 +33554432 + + + + +time_since_last_exp +me + +Titre type + + + + + + +Alpha +Ancestral +Delta +populationtrajectories + + diff --git a/tests/testthat/test-data.R b/tests/testthat/test-data.R index 57ebbaf..f1ccbd1 100644 --- a/tests/testthat/test-data.R +++ b/tests/testthat/test-data.R @@ -52,7 +52,9 @@ test_that("Natural scale data is converted to log scale for stan", { dat <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) mod <- biokinetics$new(data = dat) stan_data <- mod$get_stan_data() - expect_equal(stan_data$value, convert_log2_scale(dat, "value")$value, ignore_attr = TRUE) + expect_equal(stan_data$value, + convert_log2_scale(dat, "value")$value, + ignore_attr = TRUE) }) test_that("Log scale data is passed directly to stan", { diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R index d01a8e1..32723b3 100644 --- a/tests/testthat/test-plots.R +++ b/tests/testthat/test-plots.R @@ -8,8 +8,8 @@ test_that("Can plot prior prediction up to tmax", { 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"))] + expect_error(plot_prior_predictive(priors, data = data), "Missing required columns: time_since_last_exp") + data[, `:=`(time_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) }) @@ -28,7 +28,6 @@ test_that("Can plot prior predictions from model", { }) 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)) @@ -40,10 +39,45 @@ test_that("Prior predictions from model are the same", { vdiffr::expect_doppelganger("priorpredictive", plot) }) -test_that("Plotted data is are the same", { - skip_on_ci() +test_that("Can plot input data", { data <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) mod <- biokinetics$new(data = data) - plot <- mod$plot_data() + plot <- mod$plot_model_inputs() vdiffr::expect_doppelganger("inputdata", plot) + + mod <- biokinetics$new(data = data, covariate_formula = ~0 + infection_history) + plot <- mod$plot_model_inputs() + vdiffr::expect_doppelganger("inputdata_covariates", plot) +}) + +mock_model <- function(name, package) { + list(sample = function(x, ...) readRDS(test_path("testdata", "testdraws.rds"))) +} + +test_that("Summarised and un-summarised population trajectories give same plots", { + # note that this is using a pre-fitted model with very few iterations, so the + # fits won't look very good + local_mocked_bindings( + stan_package_model = mock_model, .package = "instantiate" + ) + mod <- biokinetics$new(file_path = system.file("delta_full.rds", package = "epikinetics"), + covariate_formula = ~0 + infection_history) + mod$fit() + trajectories <- mod$simulate_population_trajectories(summarise = TRUE) + unsummarised_trajectories <- mod$simulate_population_trajectories(summarise = FALSE) + vdiffr::expect_doppelganger("populationtrajectories", plot(trajectories)) + vdiffr::expect_doppelganger("populationtrajectories", plot(unsummarised_trajectories)) +}) + +test_that("Can plot population trajectories with data", { + local_mocked_bindings( + stan_package_model = mock_model, .package = "instantiate" + ) + # note that this is using a pre-fitted model with very few iterations, so the + # fits won't look very good + data <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) + mod <- biokinetics$new(data = data, covariate_formula = ~0 + infection_history) + fit <- mod$fit() + trajectories <- mod$simulate_population_trajectories() + vdiffr::expect_doppelganger("populationtrajectories_data", plot(trajectories, data = data)) }) From b3c8bd84e72ac6d8636f299419fade9cb50f1b1f Mon Sep 17 00:00:00 2001 From: "alex.hill@gmail.com" Date: Wed, 23 Oct 2024 23:07:58 +0100 Subject: [PATCH 04/12] fix plot test --- R/plot.R | 12 +- .../_snaps/plots/multiplecovariates.svg | 685 ++++++++++++++++++ tests/testthat/test-plots.R | 19 + 3 files changed, 712 insertions(+), 4 deletions(-) create mode 100644 tests/testthat/_snaps/plots/multiplecovariates.svg diff --git a/R/plot.R b/R/plot.R index ed11e5b..aded9bf 100644 --- a/R/plot.R +++ b/R/plot.R @@ -73,10 +73,13 @@ plot_data <- function(data, covariates) { #' a data.table containing raw data as provided to the biokinetics model. #' @export plot.biokinetics_population_trajectories <- function(x, ...) { + covariates <- attr(x, "covariates") if (!attr(x, "summarised")) { by <- setdiff(colnames(x), c("t0_pop", "tp_pop", "ts_pop", - "m1_pop", "m2_pop", "m3_pop", - ".draw", "mu")) + "m1_pop", "m2_pop", "m3_pop", + "beta_t0", "beta_tp", "beta_ts", + "beta_m1", "beta_m2", "beta_m3", + ".draw", "mu")) x <- summarise_draws( x, column_name = "mu", by = by) } @@ -101,8 +104,9 @@ plot.biokinetics_population_trajectories <- function(x, ...) { aes(x = as.integer(day - last_exp_day, units = "days"), y = value), size = 0.5, alpha = 0.5) } - plot + scale_y_continuous(trans = "log2") + - facet_wrap(eval(parse(text = facet_formula(attr(x, "covariates"))))) + + plot + + scale_y_continuous(trans = "log2") + + facet_wrap(eval(parse(text = facet_formula(covariates)))) + guides(fill = guide_legend(title = "Titre type"), colour = "none") } diff --git a/tests/testthat/_snaps/plots/multiplecovariates.svg b/tests/testthat/_snaps/plots/multiplecovariates.svg new file mode 100644 index 0000000..64d645d --- /dev/null +++ b/tests/testthat/_snaps/plots/multiplecovariates.svg @@ -0,0 +1,685 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Delta + +NA + +BNT162b2 + + + + + + + + + + +Delta + +NA + +mRNA1273 + + + + + + + + + + +Delta + +NA + +others + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Ancestral + +NA + +mRNA1273 + + + + + + + + + + +Ancestral + +NA + +others + + + + + + + + + + +Delta + +Infection naive + +NA + + + + + + + + + + +Delta + +Previously infected (Pre-Omicron) + +NA + + + + + + + + + + +Delta + +NA + +AZD1222 + + + + + + + + + + +Alpha + +NA + +others + + + + + + + + + + +Ancestral + +Infection naive + +NA + + + + + + + + + + +Ancestral + +Previously infected (Pre-Omicron) + +NA + + + + + + + + + + +Ancestral + +NA + +AZD1222 + + + + + + + + + + +Ancestral + +NA + +BNT162b2 + + + + + + + + + + +Alpha + +Infection naive + +NA + + + + + + + + + + +Alpha + +Previously infected (Pre-Omicron) + +NA + + + + + + + + + + +Alpha + +NA + +AZD1222 + + + + + + + + + + +Alpha + +NA + +BNT162b2 + + + + + + + + + + +Alpha + +NA + +mRNA1273 + + + + + + +0 +50 +100 +150 + + + + +0 +50 +100 +150 + + + + +0 +50 +100 +150 + + + + +0 +50 +100 +150 + + + + +0 +50 +100 +150 +64 +2048 +65536 +2097152 + + + + +64 +2048 +65536 +2097152 + + + + +64 +2048 +65536 +2097152 + + + + +64 +2048 +65536 +2097152 + + + + +time_since_last_exp +me + +Titre type + + + + + + +Alpha +Ancestral +Delta +multiplecovariates + + diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R index 32723b3..1074cc6 100644 --- a/tests/testthat/test-plots.R +++ b/tests/testthat/test-plots.R @@ -54,6 +54,10 @@ mock_model <- function(name, package) { list(sample = function(x, ...) readRDS(test_path("testdata", "testdraws.rds"))) } +mock_model_multiple_covariates <- function(name, package) { + list(sample = function(x, ...) readRDS(test_path("testdata", "testdraws_multiplecovariates.rds"))) +} + test_that("Summarised and un-summarised population trajectories give same plots", { # note that this is using a pre-fitted model with very few iterations, so the # fits won't look very good @@ -69,6 +73,21 @@ test_that("Summarised and un-summarised population trajectories give same plots" vdiffr::expect_doppelganger("populationtrajectories", plot(unsummarised_trajectories)) }) +test_that("Summarised and un-summarised population trajectories give same plots for multiple covariats", { + # note that this is using a pre-fitted model with very few iterations, so the + # fits won't look very good + local_mocked_bindings( + stan_package_model = mock_model_multiple_covariates, .package = "instantiate" + ) + mod <- biokinetics$new(file_path = system.file("delta_full.rds", package = "epikinetics"), + covariate_formula = ~0 + infection_history + last_vax_type) + mod$fit() + trajectories <- mod$simulate_population_trajectories(summarise = TRUE) + unsummarised_trajectories <- mod$simulate_population_trajectories(summarise = FALSE) + vdiffr::expect_doppelganger("multiplecovariates", plot(trajectories)) + vdiffr::expect_doppelganger("multiplecovariates", plot(unsummarised_trajectories)) +}) + test_that("Can plot population trajectories with data", { local_mocked_bindings( stan_package_model = mock_model, .package = "instantiate" From 76ef11a846548683b973f2f7bf5e757e8eea47c1 Mon Sep 17 00:00:00 2001 From: "alex.hill@gmail.com" Date: Wed, 23 Oct 2024 23:59:37 +0100 Subject: [PATCH 05/12] use dfault method for priors --- NAMESPACE | 2 +- R/biokinetics.R | 8 +++--- R/plot.R | 25 ++++++++++--------- ...edictive.Rd => plot.biokinetics_priors.Rd} | 10 +++++--- tests/testthat/test-plots.R | 6 ++--- 5 files changed, 27 insertions(+), 24 deletions(-) rename man/{plot_prior_predictive.Rd => plot.biokinetics_priors.Rd} (72%) diff --git a/NAMESPACE b/NAMESPACE index 7fa16cf..577ae7a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,11 +1,11 @@ # Generated by roxygen2: do not edit by hand S3method(plot,biokinetics_population_trajectories) +S3method(plot,biokinetics_priors) 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 3f6d48a..e9876f2 100644 --- a/R/biokinetics.R +++ b/R/biokinetics.R @@ -278,10 +278,10 @@ biokinetics <- R6::R6Class( #' @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) + plot(private$priors, + tmax = tmax, + n_draws = n_draws, + data = private$data) }, #' @description Plot model input data with a smoothing function. Note that #' this plot is on a log scale, regardless of whether data was provided on a diff --git a/R/plot.R b/R/plot.R index aded9bf..194a16a 100644 --- a/R/plot.R +++ b/R/plot.R @@ -4,15 +4,16 @@ #' @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 x A named list of type 'biokinetics_priors'. +#' @param \dots Further arguments passed to the method. #' @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 time_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) +plot.biokinetics_priors <- function(x, + ..., + tmax = 150, + n_draws = 2000, + data = NULL) { # Declare variables to suppress notes when compiling package # https://github.com/Rdatatable/data.table/issues/850#issuecomment-259466153 @@ -23,12 +24,12 @@ plot_prior_predictive <- function(priors, validate_required_cols(data, c("time_since_last_exp", "value")) } params <- data.table( - t0 = stats::rnorm(n_draws, priors$mu_t0, priors$sigma_t0), # titre value at t0 - tp = stats::rnorm(n_draws, priors$mu_tp, priors$sigma_tp), # time of peak - ts = stats::rnorm(n_draws, priors$mu_ts, priors$sigma_ts), # time of set point - m1 = stats::rnorm(n_draws, priors$mu_m1, priors$sigma_m1), # gradient 1 - m2 = stats::rnorm(n_draws, priors$mu_m2, priors$sigma_m2), # gradient 2 - m3 = stats::rnorm(n_draws, priors$mu_m3, priors$sigma_m3) # gradient 3 + t0 = stats::rnorm(n_draws, x$mu_t0, x$sigma_t0), # titre value at t0 + tp = stats::rnorm(n_draws, x$mu_tp, x$sigma_tp), # time of peak + ts = stats::rnorm(n_draws, x$mu_ts, x$sigma_ts), # time of set point + m1 = stats::rnorm(n_draws, x$mu_m1, x$sigma_m1), # gradient 1 + m2 = stats::rnorm(n_draws, x$mu_m2, x$sigma_m2), # gradient 2 + m3 = stats::rnorm(n_draws, x$mu_m3, x$sigma_m3) # gradient 3 ) times <- data.table(t = 1:tmax) diff --git a/man/plot_prior_predictive.Rd b/man/plot.biokinetics_priors.Rd similarity index 72% rename from man/plot_prior_predictive.Rd rename to man/plot.biokinetics_priors.Rd index fc6aa0c..3f116a8 100644 --- a/man/plot_prior_predictive.Rd +++ b/man/plot.biokinetics_priors.Rd @@ -1,14 +1,16 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/plot.R -\name{plot_prior_predictive} -\alias{plot_prior_predictive} +\name{plot.biokinetics_priors} +\alias{plot.biokinetics_priors} \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) +\method{plot}{biokinetics_priors}(x, ..., tmax = 150, n_draws = 2000, data = NULL) } \arguments{ -\item{priors}{A named list of type 'biokinetics_priors'.} +\item{x}{A named list of type 'biokinetics_priors'.} + +\item{\dots}{Further arguments passed to the method.} \item{tmax}{Integer. The number of time points in each simulated trajectory. Default 150.} diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R index 1074cc6..05fcbb2 100644 --- a/tests/testthat/test-plots.R +++ b/tests/testthat/test-plots.R @@ -1,6 +1,6 @@ test_that("Can plot prior prediction up to tmax", { priors <- biokinetics_priors() - plot <- plot_prior_predictive(priors, tmax = 100, n_draws = 500) + plot <- plot(priors, tmax = 100, n_draws = 500) expect_equal(nrow(plot$data), 100) expect_equal(length(plot$layers), 2) }) @@ -8,9 +8,9 @@ test_that("Can plot prior prediction up to tmax", { 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: time_since_last_exp") + expect_error(plot(priors, data = data), "Missing required columns: time_since_last_exp") data[, `:=`(time_since_last_exp = as.integer(day - last_exp_day, units = "days"))] - plot <- plot_prior_predictive(priors, data = data, n_draws = 500) + plot <- plot(priors, data = data, n_draws = 500) expect_equal(length(plot$layers), 3) }) From efa4635ba1c4d1b3ff11e8ebc489e8f723d0844c Mon Sep 17 00:00:00 2001 From: "alex.hill@gmail.com" Date: Thu, 24 Oct 2024 00:10:48 +0100 Subject: [PATCH 06/12] fix plot --- DESCRIPTION | 2 +- R/plot.R | 10 ++++------ man/epikinetics-package.Rd | 2 +- man/plot.biokinetics_population_trajectories.Rd | 7 ++++--- 4 files changed, 10 insertions(+), 11 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index bcab1f4..fe4cea9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -3,7 +3,7 @@ Title: Biomarker Kinetics Modelling Version: 0.0.0.9000 Authors@R: c(person("Timothy", "Russell", email = "timothy.russell@lshtm.ac.uk", role = c("aut")), person("Alex", "Hill", email = "alex.hill@gmail.com", role = c("aut", "cre"))) -Description: Fit kinetic curves to biomarker data, using a Bayesian hierarchical model +Description: Fit kinetic curves to biomarker data, using a Bayesian hierarchical model. License: GPL (>= 3) Encoding: UTF-8 Roxygen: list(markdown = TRUE) diff --git a/R/plot.R b/R/plot.R index 194a16a..0de49db 100644 --- a/R/plot.R +++ b/R/plot.R @@ -70,10 +70,10 @@ plot_data <- function(data, covariates) { #' @param x An object of class "biokinetics_population_trajectories". These are #' generated by running biokinetics$simulate_populate_trajectories(). See #' \href{../../epikinetics/html/biokinetics.html#method-biokinetics-simulate_population_trajectories}{\code{biokinetics$simulate_population_trajectories()}} -#' @param \dots Further arguments passed to the method. For now the only extra argument supported is "data", which should be -#' a data.table containing raw data as provided to the biokinetics model. +#' @param \dots Further arguments passed to the method. +#' @param data Optional data.table containing raw data as provided to the biokinetics model. #' @export -plot.biokinetics_population_trajectories <- function(x, ...) { +plot.biokinetics_population_trajectories <- function(x, ..., data = NULL) { covariates <- attr(x, "covariates") if (!attr(x, "summarised")) { by <- setdiff(colnames(x), c("t0_pop", "tp_pop", "ts_pop", @@ -84,7 +84,6 @@ plot.biokinetics_population_trajectories <- function(x, ...) { x <- summarise_draws( x, column_name = "mu", by = by) } - args <- list(...) # Declare variables to suppress notes when compiling package # https://github.com/Rdatatable/data.table/issues/850#issuecomment-259466153 time_since_last_exp <- value <- me <- titre_type <- lo <- hi <- NULL @@ -97,8 +96,7 @@ plot.biokinetics_population_trajectories <- function(x, ...) { ymax = hi, fill = titre_type), alpha = 0.5) - if (!is.null(args$data)) { - data <- args$data + if (!is.null(data)) { validate_required_cols(data) plot <- plot + geom_point(data = data, diff --git a/man/epikinetics-package.Rd b/man/epikinetics-package.Rd index 77153e1..29fb876 100644 --- a/man/epikinetics-package.Rd +++ b/man/epikinetics-package.Rd @@ -6,7 +6,7 @@ \alias{epikinetics-package} \title{epikinetics: Biomarker Kinetics Modelling} \description{ -Fit kinetic curves to biomarker data, using a Bayesian hierarchical model +Fit kinetic curves to biomarker data, using a Bayesian hierarchical model. } \author{ \strong{Maintainer}: Alex Hill \email{alex.hill@gmail.com} diff --git a/man/plot.biokinetics_population_trajectories.Rd b/man/plot.biokinetics_population_trajectories.Rd index aeecc78..6b57d13 100644 --- a/man/plot.biokinetics_population_trajectories.Rd +++ b/man/plot.biokinetics_population_trajectories.Rd @@ -4,15 +4,16 @@ \alias{plot.biokinetics_population_trajectories} \title{Plot method for "biokinetics_population_trajectories" class} \usage{ -\method{plot}{biokinetics_population_trajectories}(x, ...) +\method{plot}{biokinetics_population_trajectories}(x, ..., data = NULL) } \arguments{ \item{x}{An object of class "biokinetics_population_trajectories". These are generated by running biokinetics$simulate_populate_trajectories(). See \href{../../epikinetics/html/biokinetics.html#method-biokinetics-simulate_population_trajectories}{\code{biokinetics$simulate_population_trajectories()}}} -\item{\dots}{Further arguments passed to the method. For now the only extra argument supported is "data", which should be -a data.table containing raw data as provided to the biokinetics model.} +\item{\dots}{Further arguments passed to the method.} + +\item{data}{Optional data.table containing raw data as provided to the biokinetics model.} } \description{ Plot method for "biokinetics_population_trajectories" class From 088629c7b315edd5f1092a5aec0f6c6604549a30 Mon Sep 17 00:00:00 2001 From: "alex.hill@gmail.com" Date: Thu, 24 Oct 2024 09:39:30 +0100 Subject: [PATCH 07/12] don't use quantiles for unsummarised --- R/biokinetics.R | 15 +- R/plot.R | 27 +- .../_snaps/plots/multiplecovariates-unsum.svg | 620 ++++ .../_snaps/plots/multiplecovariates.svg | 66 +- .../plots/populationtrajectories-data.svg | 2670 ++++++++--------- .../plots/populationtrajectories-unsum.svg | 216 ++ .../_snaps/plots/populationtrajectories.svg | 72 +- tests/testthat/test-plots.R | 14 +- 8 files changed, 2270 insertions(+), 1430 deletions(-) create mode 100644 tests/testthat/_snaps/plots/multiplecovariates-unsum.svg create mode 100644 tests/testthat/_snaps/plots/populationtrajectories-unsum.svg diff --git a/R/biokinetics.R b/R/biokinetics.R index e9876f2..fe9e1e6 100644 --- a/R/biokinetics.R +++ b/R/biokinetics.R @@ -142,11 +142,6 @@ biokinetics <- R6::R6Class( dt_out[, t_id := NULL] - if (private$scale == "natural") { - dt_out <- convert_log2_scale_inverse( - dt_out, vars_to_transform = "mu") - } - if (summarise == TRUE) { logger::log_info("Summarising into quantiles") dt_out <- summarise_draws( @@ -414,6 +409,16 @@ biokinetics <- R6::R6Class( dt_out <- dt_out[ , lapply(.SD, function(x) if (is.factor(x)) forcats::fct_drop(x) else x)] + if (private$scale == "natural") { + if (summarise) { + dt_out <- convert_log2_scale_inverse( + dt_out, vars_to_transform = c("me", "lo", "hi")) + } else { + dt_out <- convert_log2_scale_inverse( + dt_out, vars_to_transform = "mu") + } + } + class(dt_out) <- append("biokinetics_population_trajectories", class(dt_out)) attr(dt_out, "summarised") <- summarise attr(dt_out, "covariates") <- private$all_formula_vars diff --git a/R/plot.R b/R/plot.R index 0de49db..0d82e15 100644 --- a/R/plot.R +++ b/R/plot.R @@ -75,26 +75,23 @@ plot_data <- function(data, covariates) { #' @export plot.biokinetics_population_trajectories <- function(x, ..., data = NULL) { covariates <- attr(x, "covariates") - if (!attr(x, "summarised")) { - by <- setdiff(colnames(x), c("t0_pop", "tp_pop", "ts_pop", - "m1_pop", "m2_pop", "m3_pop", - "beta_t0", "beta_tp", "beta_ts", - "beta_m1", "beta_m2", "beta_m3", - ".draw", "mu")) - x <- summarise_draws( - x, column_name = "mu", by = by) - } + # Declare variables to suppress notes when compiling package # https://github.com/Rdatatable/data.table/issues/850#issuecomment-259466153 time_since_last_exp <- value <- me <- titre_type <- lo <- hi <- NULL day <- last_exp_day <- NULL - plot <- ggplot(x) + - geom_line(aes(x = time_since_last_exp, y = me, colour = titre_type)) + - geom_ribbon(aes(x = time_since_last_exp, - ymin = lo, - ymax = hi, - fill = titre_type), alpha = 0.5) + if (attr(x, "summarised")) { + plot <- ggplot(x) + + geom_line(aes(x = time_since_last_exp, y = me, colour = titre_type)) + + geom_ribbon(aes(x = time_since_last_exp, + ymin = lo, + ymax = hi, + fill = titre_type), alpha = 0.5) + } else { + plot <- ggplot(x) + + geom_line(aes(x = time_since_last_exp, y = mu, colour = titre_type), alpha = 0.5) + } if (!is.null(data)) { validate_required_cols(data) diff --git a/tests/testthat/_snaps/plots/multiplecovariates-unsum.svg b/tests/testthat/_snaps/plots/multiplecovariates-unsum.svg new file mode 100644 index 0000000..63f4a71 --- /dev/null +++ b/tests/testthat/_snaps/plots/multiplecovariates-unsum.svg @@ -0,0 +1,620 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Delta + +NA + +BNT162b2 + + + + + + + + + + +Delta + +NA + +mRNA1273 + + + + + + + + + + +Delta + +NA + +others + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Ancestral + +NA + +mRNA1273 + + + + + + + + + + +Ancestral + +NA + +others + + + + + + + + + + +Delta + +Infection naive + +NA + + + + + + + + + + +Delta + +Previously infected (Pre-Omicron) + +NA + + + + + + + + + + +Delta + +NA + +AZD1222 + + + + + + + + + + +Alpha + +NA + +others + + + + + + + + + + +Ancestral + +Infection naive + +NA + + + + + + + + + + +Ancestral + +Previously infected (Pre-Omicron) + +NA + + + + + + + + + + +Ancestral + +NA + +AZD1222 + + + + + + + + + + +Ancestral + +NA + +BNT162b2 + + + + + + + + + + +Alpha + +Infection naive + +NA + + + + + + + + + + +Alpha + +Previously infected (Pre-Omicron) + +NA + + + + + + + + + + +Alpha + +NA + +AZD1222 + + + + + + + + + + +Alpha + +NA + +BNT162b2 + + + + + + + + + + +Alpha + +NA + +mRNA1273 + + + + + + +0 +50 +100 +150 + + + + +0 +50 +100 +150 + + + + +0 +50 +100 +150 + + + + +0 +50 +100 +150 + + + + +0 +50 +100 +150 +64 +2048 +65536 +2097152 + + + + +64 +2048 +65536 +2097152 + + + + +64 +2048 +65536 +2097152 + + + + +64 +2048 +65536 +2097152 + + + + +time_since_last_exp +mu +multiplecovariates_unsum + + diff --git a/tests/testthat/_snaps/plots/multiplecovariates.svg b/tests/testthat/_snaps/plots/multiplecovariates.svg index 64d645d..61c4096 100644 --- a/tests/testthat/_snaps/plots/multiplecovariates.svg +++ b/tests/testthat/_snaps/plots/multiplecovariates.svg @@ -27,7 +27,7 @@ - + @@ -43,9 +43,9 @@ - + - + @@ -57,7 +57,7 @@ - + @@ -72,7 +72,7 @@ - + @@ -87,7 +87,7 @@ - + @@ -102,9 +102,9 @@ - - - + + + @@ -117,10 +117,10 @@ - - + + - + @@ -132,7 +132,7 @@ - + @@ -147,10 +147,10 @@ - - - - + + + + @@ -162,7 +162,7 @@ - + @@ -177,7 +177,7 @@ - + @@ -192,7 +192,7 @@ - + @@ -207,10 +207,10 @@ - - + + - + @@ -222,9 +222,9 @@ - - - + + + @@ -237,7 +237,7 @@ - + @@ -262,8 +262,8 @@ - - + + @@ -276,9 +276,9 @@ - - - + + + @@ -291,7 +291,7 @@ - + diff --git a/tests/testthat/_snaps/plots/populationtrajectories-data.svg b/tests/testthat/_snaps/plots/populationtrajectories-data.svg index 6340bda..5b5814f 100644 --- a/tests/testthat/_snaps/plots/populationtrajectories-data.svg +++ b/tests/testthat/_snaps/plots/populationtrajectories-data.svg @@ -27,409 +27,409 @@ - - - - + + + + - + - + - - - + + + - + - - - - - - - - + + + + + + + + - - - - + + + + - + - - + + - + - + - + - - + + - - - + + + - - + + - + - - - - - - - - + + + + + + + + - + - + - - + + - - - + + + - - + + - - - - + + + + - - - + + + - - + + - + - - - - - + + + + + - + - - + + - + - + - + - + - - + + - + - - - - + + + + - + - - + + - + - - - + + + - + - - + + - - - - - + + + + + - + - + - + - - - - - - + + + + + + - - - - - + + + + + - + - + - - + + - + - - + + - + - + - + - - + + - - - - + + + + - - - + + + - - - - - + + + + + - + - - - - - + + + + + - - + + - - + + - + - + - - - - + + + + - - - - - - - - - - - - + + + + + + + + + + + + - - - + + + - - + + - - - - - - - + + + + + + + - + - - - - - + + + + + - - + + - - + + - - - - + + + + - - + + - + - - + + - - - + + + - - - - - + + + + + - - - + + + - - + + - + - + - - + + - + - - - - - + + + + + - + - + - - - + + + - + - - - + + + - - + + - - + + - - + + - - - - - + + + + + @@ -437,149 +437,149 @@ - + - - - - - + + + + + - + - + - - - - - - - + + + + + + + - + - - + + - - + + - - - + + + - + - + - + - - - + + + - - + + - + - + - - - - + + + + - - - + + + - - - - - + + + + + - - + + - + - - - - + + + + - - - - + + + + - - + + - - - - - + + + + + - + - - - - + + + + - + - + - - + + - + - - - - - - + + + + + + - - - + + + - + - + - + - + @@ -591,209 +591,209 @@ - - - + + + - - - + + + - - - + + + - - - - + + + + - + - - + + - - - + + + - + - - - - - - - + + + + + + + - - - - + + + + - - - - - - + + + + + + - + - - - - - - - + + + + + + + - - - + + + - + - - - + + + - + - - + + - + - + - - + + - - + + - + - - - - - - - - - - + + + + + + + + + + - - - - + + + + - + - - - - - - - - - - + + + + + + + + + + - - + + - - - + + + - + - - - - + + + + - + - - - - - - - - + + + + + + + + - - + + - - - + + + - - - - - - - - + + + + + + + + - - - + + + - - - - - - + + + + + + - + - - - - + + + + - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + @@ -805,214 +805,214 @@ - - - - - - - + + + + + + + - - - + + + - - - - + + + + - - - + + + - + - - - + + + - - + + - + - - - - + + + + - - - - - - + + + + + + - - - - - - + + + + + + - - + + - - - - + + + + - + - + - + - - - - - - - - - + + + + + + + + + - - - - - - + + + + + + - + - - - + + + - - - - - - + + + + + + - - + + - - + + - + - - - + + + - + - + - - + + - - + + - + - - + + - + - + - - - - + + + + - + - - + + - + - - - - - - - - - - + + + + + + + + + + - + - - - + + + - + - + - - - + + + - - + + - - - + + + - - - - - - - + + + + + + + - - - + + + @@ -1025,62 +1025,62 @@ - - - - + + + + - + - + - + - + - + - + - - + + - + - - + + - - - - + + + + - + - + - + - + @@ -1088,65 +1088,65 @@ - + - - + + - - - - + + + + - + - + - + - - + + - + - - - - - - + + + + + + - - + + - + - + - + - - - - - + + + + + - + - - + + - + - + @@ -1155,304 +1155,304 @@ - - + + - + - + - - - - + + + + - + - + - - - - - + + + + + - - + + - - - - + + + + - - - + + + - + - + - - + + - + - - - - + + + + - - - + + + - - + + - + - + - + - + - - + + - + - + - - + + - - - + + + - + - - + + - + - - + + - + - - - + + + - - - - + + + + - + - + - + - - - - - + + + + + - + - - + + - + - + - + - - - + + + - + - + - + - - - + + + - - - - - - + + + + + + - - + + - - - - + + + + - - - - - + + + + + - + - - + + - + - + - - + + - + - - - + + + - - + + - + - - - - + + + + - + - - + + - - - - - - + + + + + + - - - + + + - + - + - + - + - - + + - - - + + + - - - + + + - - - + + + - + - + - + @@ -1463,124 +1463,124 @@ - - + + - - - - + + + + - + - + - + - - + + - + - - - + + + - + - - + + - - - - - + + + + + - + - - - + + + - + - - + + - + - + - - - - - + + + + + - - + + - - - - + + + + - + - - - + + + - + - - + + - - + + - + - + @@ -1596,399 +1596,399 @@ - - - - - + + + + + - - + + - - - + + + - - - - - - + + + + + + - + - - - - - + + + + + - + - - - - - + + + + + - - - - + + + + - + - - - + + + - - + + - + - - - + + + - - - + + + - - - + + + - - - - + + + + - - - + + + - + - - - - - + + + + + - + - + - - - - + + + + - + - + - + - + - - + + - - - + + + - + - - - + + + - - - + + + - + - + - - - + + + - - - - + + + + - - - - + + + + - - - - + + + + - - - - - + + + + + - - - - + + + + - + - - - - - - + + + + + + - - + + - + - + - - + + - - - - - + + + + + - - - - + + + + - + - - + + - + - + - + - - + + - - + + - + - + - + - - - - + + + + - - + + - - - + + + - - + + - - + + - + - - - - - + + + + + - - - - + + + + - - + + - - - - + + + + - - - + + + - - - + + + - + - - - - - - - - + + + + + + + + - + - - + + - - - + + + - - - - - + + + + + - + - + - - - - + + + + - - + + - - - + + + - + - - - - - + + + + + - + - + - + - - - - + + + + - - - + + + - - - - - - + + + + + + - - + + - - - - + + + + - + @@ -1998,140 +1998,140 @@ - - - - + + + + - - + + - - + + - - + + - + - - - - + + + + - + - + - - - - - - - - + + + + + + + + - + - + - - + + - + - - - - + + + + - + - + - - - - - + + + + + - + - - + + - - - - - + + + + + - + - - - - + + + + - + - + - + - - + + - - - + + + - - + + - - - + + + - + - - - + + + - + - + - + @@ -2146,221 +2146,221 @@ - - - - - - + + + + + + - + - - - - + + + + - - + + - + - - + + - - - + + + - + - - - + + + - + - - + + - + - - - - - + + + + + - + - - - - - + + + + + - - - - + + + + - + - - - + + + - - + + - - + + - - + + - + - + - + - - - - - - - - - + + + + + + + + + - - - + + + - - + + - - - - - + + + + + - - + + - - - + + + - + - - + + - + - + - - - - - - + + + + + + - + - - - - - + + + + + - - + + - - - - + + + + - + - - + + - + - + - - - - - + + + + + - - - - + + + + - + - - - - - - - - + + + + + + + + - - - - - - - - - + + + + + + + + + - + @@ -2466,22 +2466,22 @@ 200 400 600 -128 -8192 -524288 -33554432 +128 +8192 +524288 +33554432 - - - + + + 128 -8192 -524288 -33554432 +8192 +524288 +33554432 - - - + + + time_since_last_exp me diff --git a/tests/testthat/_snaps/plots/populationtrajectories-unsum.svg b/tests/testthat/_snaps/plots/populationtrajectories-unsum.svg new file mode 100644 index 0000000..43579e1 --- /dev/null +++ b/tests/testthat/_snaps/plots/populationtrajectories-unsum.svg @@ -0,0 +1,216 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Ancestral + +Previously infected (Pre-Omicron) + + + + + + + + + + +Delta + +Infection naive + + + + + + + + + + +Delta + +Previously infected (Pre-Omicron) + + + + + + + + + + +Alpha + +Infection naive + + + + + + + + + + +Alpha + +Previously infected (Pre-Omicron) + + + + + + + + + + +Ancestral + +Infection naive + + + + + + +0 +50 +100 +150 + + + + +0 +50 +100 +150 + + + + +0 +50 +100 +150 +128 +8192 +524288 +33554432 + + + + +128 +8192 +524288 +33554432 + + + + +time_since_last_exp +mu +populationtrajectories_unsum + + diff --git a/tests/testthat/_snaps/plots/populationtrajectories.svg b/tests/testthat/_snaps/plots/populationtrajectories.svg index 3c2d3b9..1aa7f3f 100644 --- a/tests/testthat/_snaps/plots/populationtrajectories.svg +++ b/tests/testthat/_snaps/plots/populationtrajectories.svg @@ -27,10 +27,10 @@ - - - - + + + + @@ -42,9 +42,9 @@ - - - + + + @@ -57,10 +57,10 @@ - - - - + + + + @@ -72,10 +72,10 @@ - - - - + + + + @@ -87,10 +87,10 @@ - - - - + + + + @@ -102,10 +102,10 @@ - - - - + + + + @@ -211,22 +211,22 @@ 50 100 150 -128 -8192 -524288 -33554432 +128 +8192 +524288 +33554432 - - - + + + 128 -8192 -524288 -33554432 +8192 +524288 +33554432 - - - + + + time_since_last_exp me diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R index 05fcbb2..5746122 100644 --- a/tests/testthat/test-plots.R +++ b/tests/testthat/test-plots.R @@ -58,7 +58,7 @@ mock_model_multiple_covariates <- function(name, package) { list(sample = function(x, ...) readRDS(test_path("testdata", "testdraws_multiplecovariates.rds"))) } -test_that("Summarised and un-summarised population trajectories give same plots", { +test_that("Cna plot summarised and un-summarised population trajectories", { # note that this is using a pre-fitted model with very few iterations, so the # fits won't look very good local_mocked_bindings( @@ -68,12 +68,13 @@ test_that("Summarised and un-summarised population trajectories give same plots" covariate_formula = ~0 + infection_history) mod$fit() trajectories <- mod$simulate_population_trajectories(summarise = TRUE) - unsummarised_trajectories <- mod$simulate_population_trajectories(summarise = FALSE) vdiffr::expect_doppelganger("populationtrajectories", plot(trajectories)) - vdiffr::expect_doppelganger("populationtrajectories", plot(unsummarised_trajectories)) + + unsummarised_trajectories <- mod$simulate_population_trajectories(summarise = FALSE) + vdiffr::expect_doppelganger("populationtrajectories_unsum", plot(unsummarised_trajectories)) }) -test_that("Summarised and un-summarised population trajectories give same plots for multiple covariats", { +test_that("Can plot summarised and un-summarised population trajectories for multiple covariates", { # note that this is using a pre-fitted model with very few iterations, so the # fits won't look very good local_mocked_bindings( @@ -83,9 +84,10 @@ test_that("Summarised and un-summarised population trajectories give same plots covariate_formula = ~0 + infection_history + last_vax_type) mod$fit() trajectories <- mod$simulate_population_trajectories(summarise = TRUE) - unsummarised_trajectories <- mod$simulate_population_trajectories(summarise = FALSE) vdiffr::expect_doppelganger("multiplecovariates", plot(trajectories)) - vdiffr::expect_doppelganger("multiplecovariates", plot(unsummarised_trajectories)) + + unsummarised_trajectories <- mod$simulate_population_trajectories(summarise = FALSE) + vdiffr::expect_doppelganger("multiplecovariates_unsum", plot(unsummarised_trajectories)) }) test_that("Can plot population trajectories with data", { From 9c40c2d55779099455fc6ff447ee3477f7e5a2d6 Mon Sep 17 00:00:00 2001 From: "alex.hill@gmail.com" Date: Tue, 29 Oct 2024 10:13:53 +0000 Subject: [PATCH 08/12] pr suggestions --- R/biokinetics.R | 2 +- R/plot.R | 7 +- .../_snaps/plots/multiplecovariates-unsum.svg | 2898 ++++++++++++++++- .../plots/populationtrajectories-unsum.svg | 966 +++++- 4 files changed, 3845 insertions(+), 28 deletions(-) diff --git a/R/biokinetics.R b/R/biokinetics.R index fe9e1e6..8cceffc 100644 --- a/R/biokinetics.R +++ b/R/biokinetics.R @@ -283,7 +283,7 @@ biokinetics <- R6::R6Class( #' log or a natural scale. #' @return A ggplot2 object. plot_model_inputs = function() { - plot_data(private$data, private$all_formula_vars) + plot_sero_data(private$data, private$all_formula_vars) }, #' @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. diff --git a/R/plot.R b/R/plot.R index 0d82e15..4890f64 100644 --- a/R/plot.R +++ b/R/plot.R @@ -52,7 +52,7 @@ plot.biokinetics_priors <- function(x, plot } -plot_data <- function(data, covariates) { +plot_sero_data <- function(data, covariates = character(0)) { validate_required_cols(data, c("time_since_last_exp", "value", "titre_type")) # Declare variables to suppress notes when compiling package # https://github.com/Rdatatable/data.table/issues/850#issuecomment-259466153 @@ -68,7 +68,7 @@ plot_data <- function(data, covariates) { #' Plot method for "biokinetics_population_trajectories" class #' #' @param x An object of class "biokinetics_population_trajectories". These are -#' generated by running biokinetics$simulate_populate_trajectories(). See +#' generated by running biokinetics$simulate_population_trajectories(). See #' \href{../../epikinetics/html/biokinetics.html#method-biokinetics-simulate_population_trajectories}{\code{biokinetics$simulate_population_trajectories()}} #' @param \dots Further arguments passed to the method. #' @param data Optional data.table containing raw data as provided to the biokinetics model. @@ -90,7 +90,8 @@ plot.biokinetics_population_trajectories <- function(x, ..., data = NULL) { fill = titre_type), alpha = 0.5) } else { plot <- ggplot(x) + - geom_line(aes(x = time_since_last_exp, y = mu, colour = titre_type), alpha = 0.5) + geom_line(aes(x = time_since_last_exp, y = mu, + colour = titre_type, group = .draw), alpha = 0.1, linewidth = 0.1) } if (!is.null(data)) { diff --git a/tests/testthat/_snaps/plots/multiplecovariates-unsum.svg b/tests/testthat/_snaps/plots/multiplecovariates-unsum.svg index 63f4a71..b459444 100644 --- a/tests/testthat/_snaps/plots/multiplecovariates-unsum.svg +++ b/tests/testthat/_snaps/plots/multiplecovariates-unsum.svg @@ -27,7 +27,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -39,7 +198,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -51,7 +369,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -63,7 +540,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -75,7 +711,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -87,7 +882,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -99,7 +1053,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -111,7 +1224,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -123,7 +1395,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -135,7 +1566,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -147,7 +1737,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -159,7 +1908,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -171,7 +2079,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -183,7 +2250,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -195,7 +2421,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -216,7 +2601,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -228,7 +2772,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -240,7 +2943,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/_snaps/plots/populationtrajectories-unsum.svg b/tests/testthat/_snaps/plots/populationtrajectories-unsum.svg index 43579e1..f06ec39 100644 --- a/tests/testthat/_snaps/plots/populationtrajectories-unsum.svg +++ b/tests/testthat/_snaps/plots/populationtrajectories-unsum.svg @@ -27,7 +27,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -39,7 +198,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -51,7 +369,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -63,7 +540,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -75,7 +711,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -87,7 +882,166 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 61ac60ea9c39cae8e2313488fba306d2675c5c49 Mon Sep 17 00:00:00 2001 From: "alex.hill@gmail.com" Date: Tue, 29 Oct 2024 10:16:17 +0000 Subject: [PATCH 09/12] make n_draws consistent --- R/biokinetics.R | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/R/biokinetics.R b/R/biokinetics.R index 8cceffc..980493b 100644 --- a/R/biokinetics.R +++ b/R/biokinetics.R @@ -191,7 +191,7 @@ biokinetics <- R6::R6Class( return(dt) }, - extract_parameters = function(params, n_draws = 2500) { + extract_parameters = function(params, n_draws) { private$check_fitted() params_proc <- rlang::parse_exprs(params) @@ -306,9 +306,9 @@ biokinetics <- R6::R6Class( }, #' @description Extract fitted population parameters #' @return A data.table - #' @param n_draws Numeric + #' @param n_draws Integer. Default 2000. #' @param human_readable_covariates Logical. Default TRUE. - extract_population_parameters = function(n_draws = 2500, + extract_population_parameters = function(n_draws = 2000, human_readable_covariates = TRUE) { private$check_fitted() has_covariates <- length(private$all_formula_vars) > 0 @@ -342,10 +342,10 @@ biokinetics <- R6::R6Class( }, #' @description Extract fitted individual parameters #' @return A data.table - #' @param n_draws Numeric + #' @param n_draws Integer. Default 2000. #' @param include_variation_params Logical #' @param human_readable_covariates Logical. Default TRUE. - extract_individual_parameters = function(n_draws = 2500, + extract_individual_parameters = function(n_draws = 2000, include_variation_params = TRUE, human_readable_covariates = TRUE) { private$check_fitted() @@ -387,11 +387,11 @@ biokinetics <- R6::R6Class( #' @param summarise Boolean. Default TRUE. If TRUE, summarises over draws from posterior parameter distributions to #' return 0.025, 0.5 and 0.975 quantiles, labelled lo, me and hi, respectively. If FALSE returns values for individual #' draws from posterior parameter distributions. - #' @param n_draws Integer. Maximum number of samples to include. Default 2500. + #' @param n_draws Integer. Maximum number of samples to include. Default 2000. simulate_population_trajectories = function( t_max = 150, summarise = TRUE, - n_draws = 2500) { + n_draws = 2000) { private$check_fitted() validate_numeric(t_max) @@ -428,8 +428,8 @@ biokinetics <- R6::R6Class( #' @return A data.table of peak and set titre values. Columns are tire_type, mu_p, mu_s, rel_drop_me, mu_p_me, #' mu_s_me, and a column for each covariate. See the data vignette for details: #' \code{vignette("data", package = "epikinetics")} - #' @param n_draws Integer. Maximum number of samples to include. Default 2500. - population_stationary_points = function(n_draws = 2500) { + #' @param n_draws Integer. Maximum number of samples to include. Default 2000. + population_stationary_points = function(n_draws = 2000) { private$check_fitted() validate_numeric(n_draws) @@ -481,11 +481,11 @@ biokinetics <- R6::R6Class( #' @param summarise Boolean. If TRUE, average the individual trajectories to get lo, me and #' hi values for the population, disaggregated by titre type. If FALSE return the indidivudal trajectories. #' Default TRUE. - #' @param n_draws Integer. Maximum number of samples to draw. Default 2500. + #' @param n_draws Integer. Maximum number of samples to draw. Default 2000. #' @param time_shift Integer. Number of days to adjust the exposure day by. Default 0. simulate_individual_trajectories = function( summarise = TRUE, - n_draws = 2500, + n_draws = 2000, time_shift = 0) { private$check_fitted() validate_logical(summarise) From cd21e6e003e07a59a28638dfea4cd815e82af5c5 Mon Sep 17 00:00:00 2001 From: "alex.hill@gmail.com" Date: Tue, 29 Oct 2024 11:28:26 +0000 Subject: [PATCH 10/12] update docs --- DESCRIPTION | 2 +- man/biokinetics.Rd | 20 +++++++++---------- ...lot.biokinetics_population_trajectories.Rd | 2 +- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index fe4cea9..81b40d8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -7,7 +7,7 @@ Description: Fit kinetic curves to biomarker data, using a Bayesian hierarchical License: GPL (>= 3) Encoding: UTF-8 Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 Imports: cmdstanr, data.table, diff --git a/man/biokinetics.Rd b/man/biokinetics.Rd index feb9501..4e5f3f0 100644 --- a/man/biokinetics.Rd +++ b/man/biokinetics.Rd @@ -156,7 +156,7 @@ A CmdStanMCMC fitted model object: \url{https://mc-stan.org/cmdstanr/reference/C Extract fitted population parameters \subsection{Usage}{ \if{html}{\out{
}}\preformatted{biokinetics$extract_population_parameters( - n_draws = 2500, + n_draws = 2000, human_readable_covariates = TRUE )}\if{html}{\out{
}} } @@ -164,7 +164,7 @@ Extract fitted population parameters \subsection{Arguments}{ \if{html}{\out{
}} \describe{ -\item{\code{n_draws}}{Numeric} +\item{\code{n_draws}}{Integer. Default 2000.} \item{\code{human_readable_covariates}}{Logical. Default TRUE.} } @@ -181,7 +181,7 @@ A data.table Extract fitted individual parameters \subsection{Usage}{ \if{html}{\out{
}}\preformatted{biokinetics$extract_individual_parameters( - n_draws = 2500, + n_draws = 2000, include_variation_params = TRUE, human_readable_covariates = TRUE )}\if{html}{\out{
}} @@ -190,7 +190,7 @@ Extract fitted individual parameters \subsection{Arguments}{ \if{html}{\out{
}} \describe{ -\item{\code{n_draws}}{Numeric} +\item{\code{n_draws}}{Integer. Default 2000.} \item{\code{include_variation_params}}{Logical} @@ -211,7 +211,7 @@ Process the model results into a data table of titre values over time. \if{html}{\out{
}}\preformatted{biokinetics$simulate_population_trajectories( t_max = 150, summarise = TRUE, - n_draws = 2500 + n_draws = 2000 )}\if{html}{\out{
}} } @@ -224,7 +224,7 @@ Process the model results into a data table of titre values over time. return 0.025, 0.5 and 0.975 quantiles, labelled lo, me and hi, respectively. If FALSE returns values for individual draws from posterior parameter distributions.} -\item{\code{n_draws}}{Integer. Maximum number of samples to include. Default 2500.} +\item{\code{n_draws}}{Integer. Maximum number of samples to include. Default 2000.} } \if{html}{\out{
}} } @@ -242,13 +242,13 @@ titre_type and a column for each covariate in the hierarchical model. See the da \subsection{Method \code{population_stationary_points()}}{ Process the stan model results into a data.table. \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{biokinetics$population_stationary_points(n_draws = 2500)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{biokinetics$population_stationary_points(n_draws = 2000)}\if{html}{\out{
}} } \subsection{Arguments}{ \if{html}{\out{
}} \describe{ -\item{\code{n_draws}}{Integer. Maximum number of samples to include. Default 2500.} +\item{\code{n_draws}}{Integer. Maximum number of samples to include. Default 2000.} } \if{html}{\out{
}} } @@ -267,7 +267,7 @@ computationally expensive and may take a while to run if n_draws is large. \subsection{Usage}{ \if{html}{\out{
}}\preformatted{biokinetics$simulate_individual_trajectories( summarise = TRUE, - n_draws = 2500, + n_draws = 2000, time_shift = 0 )}\if{html}{\out{
}} } @@ -279,7 +279,7 @@ computationally expensive and may take a while to run if n_draws is large. hi values for the population, disaggregated by titre type. If FALSE return the indidivudal trajectories. Default TRUE.} -\item{\code{n_draws}}{Integer. Maximum number of samples to draw. Default 2500.} +\item{\code{n_draws}}{Integer. Maximum number of samples to draw. Default 2000.} \item{\code{time_shift}}{Integer. Number of days to adjust the exposure day by. Default 0.} } diff --git a/man/plot.biokinetics_population_trajectories.Rd b/man/plot.biokinetics_population_trajectories.Rd index 6b57d13..333ae60 100644 --- a/man/plot.biokinetics_population_trajectories.Rd +++ b/man/plot.biokinetics_population_trajectories.Rd @@ -8,7 +8,7 @@ } \arguments{ \item{x}{An object of class "biokinetics_population_trajectories". These are -generated by running biokinetics$simulate_populate_trajectories(). See +generated by running biokinetics$simulate_population_trajectories(). See \href{../../epikinetics/html/biokinetics.html#method-biokinetics-simulate_population_trajectories}{\code{biokinetics$simulate_population_trajectories()}}} \item{\dots}{Further arguments passed to the method.} From ea4ed856d86f724b273e18f005c7dba30ddd4826 Mon Sep 17 00:00:00 2001 From: "alex.hill@gmail.com" Date: Wed, 30 Oct 2024 15:15:28 +0000 Subject: [PATCH 11/12] condiitonal log transform and shape based on censoring --- R/biokinetics.R | 1 + R/plot.R | 29 +- R/utils.R | 1 - src/stan/antibody_kinetics_main.stan | 2 +- .../plots/populationtrajectories-data.svg | 4540 ++++++++-------- .../plots/populationtrajectories-logscale.svg | 2512 +++++++++ .../testthat/_snaps/plots/priorpredictive.svg | 4551 +++++++++-------- tests/testthat/test-plots.R | 22 +- 8 files changed, 7111 insertions(+), 4547 deletions(-) create mode 100644 tests/testthat/_snaps/plots/populationtrajectories-logscale.svg diff --git a/R/biokinetics.R b/R/biokinetics.R index 980493b..e93e1e1 100644 --- a/R/biokinetics.R +++ b/R/biokinetics.R @@ -421,6 +421,7 @@ biokinetics <- R6::R6Class( class(dt_out) <- append("biokinetics_population_trajectories", class(dt_out)) attr(dt_out, "summarised") <- summarise + attr(dt_out, "scale") <- private$scale attr(dt_out, "covariates") <- private$all_formula_vars dt_out }, diff --git a/R/plot.R b/R/plot.R index 4890f64..6cc2a33 100644 --- a/R/plot.R +++ b/R/plot.R @@ -47,7 +47,18 @@ plot.biokinetics_priors <- function(x, geom_ribbon(aes(x = t, ymin = lo, ymax = hi), alpha = 0.5) if (!is.null(data)) { - plot <- plot + geom_point(data = data, aes(x = time_since_last_exp, y = value)) + if (!("censored" %in% colnames(data))) { + data$censored <- FALSE + } else { + data$censored <- data$censored != 0 + } + dat <- data[time_since_last_exp <= tmax,] + plot <- plot + + geom_point(data = dat, size = 0.5, + aes(x = time_since_last_exp, + y = value, + shape = censored)) + + guides(shape = guide_legend(title = "Censored")) } plot } @@ -93,18 +104,26 @@ plot.biokinetics_population_trajectories <- function(x, ..., data = NULL) { geom_line(aes(x = time_since_last_exp, y = mu, colour = titre_type, group = .draw), alpha = 0.1, linewidth = 0.1) } - if (!is.null(data)) { validate_required_cols(data) + if (!("censored" %in% colnames(data))) { + data$censored <- FALSE + } else { + data$censored <- data$censored != 0 + } plot <- plot + geom_point(data = data, aes(x = as.integer(day - last_exp_day, units = "days"), - y = value), size = 0.5, alpha = 0.5) + y = value, shape = censored), size = 0.5, alpha = 0.5) + } + if (attr(x, "scale") == "natural") { + plot <- plot + scale_y_continuous(trans = "log2") } plot + - scale_y_continuous(trans = "log2") + facet_wrap(eval(parse(text = facet_formula(covariates)))) + - guides(fill = guide_legend(title = "Titre type"), colour = "none") + guides(fill = guide_legend(title = "Titre type"), + colour = "none", + shape = guide_legend(title = "Censored")) } facet_formula <- function(covariates) { diff --git a/R/utils.R b/R/utils.R index 6d770f4..532f9c8 100644 --- a/R/utils.R +++ b/R/utils.R @@ -1,7 +1,6 @@ convert_log2_scale <- function( dt_in, vars_to_transform = "value", simplify_limits = TRUE) { - dt_out <- data.table::copy(dt_in) for (var in vars_to_transform) { if (simplify_limits == TRUE) { diff --git a/src/stan/antibody_kinetics_main.stan b/src/stan/antibody_kinetics_main.stan index 808fc88..61708a2 100644 --- a/src/stan/antibody_kinetics_main.stan +++ b/src/stan/antibody_kinetics_main.stan @@ -149,7 +149,7 @@ model { // Likelihood for observations at lower limit: log2(value/5) = 0 target += normal_lcdf(0 | mu[cens_lo_idx], sigma); - // Censoring at log2(value/5) = 7, originally value = 2560 + // Censoring at log2(value/5) = 9, originally value = 2560 target += normal_lccdf(9 | mu[cens_hi_idx], sigma); // Covariate-level mean priors, parameterised from previous studies diff --git a/tests/testthat/_snaps/plots/populationtrajectories-data.svg b/tests/testthat/_snaps/plots/populationtrajectories-data.svg index 5b5814f..25a0e04 100644 --- a/tests/testthat/_snaps/plots/populationtrajectories-data.svg +++ b/tests/testthat/_snaps/plots/populationtrajectories-data.svg @@ -31,555 +31,555 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -595,205 +595,205 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -809,211 +809,211 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -1029,562 +1029,562 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -1600,541 +1600,541 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -2150,217 +2150,217 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -2484,17 +2484,25 @@ time_since_last_exp me - -Titre type - - - - - - -Alpha -Ancestral -Delta + +Titre type + + + + + + +Alpha +Ancestral +Delta + +Censored + + + + +FALSE +TRUE populationtrajectories_data diff --git a/tests/testthat/_snaps/plots/populationtrajectories-logscale.svg b/tests/testthat/_snaps/plots/populationtrajectories-logscale.svg new file mode 100644 index 0000000..dceb8df --- /dev/null +++ b/tests/testthat/_snaps/plots/populationtrajectories-logscale.svg @@ -0,0 +1,2512 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Ancestral + +Previously infected (Pre-Omicron) + + + + + + + + + + +Delta + +Infection naive + + + + + + + + + + +Delta + +Previously infected (Pre-Omicron) + + + + + + + + + + +Alpha + +Infection naive + + + + + + + + + + +Alpha + +Previously infected (Pre-Omicron) + + + + + + + + + + +Ancestral + +Infection naive + + + + + + +0 +200 +400 +600 + + + + +0 +200 +400 +600 + + + + +0 +200 +400 +600 +0 +5 +10 +15 +20 + + + + + +0 +5 +10 +15 +20 + + + + + +time_since_last_exp +me + +Titre type + + + + + + +Alpha +Ancestral +Delta + +Censored + + + + +FALSE +TRUE +populationtrajectories_logscale + + diff --git a/tests/testthat/_snaps/plots/priorpredictive.svg b/tests/testthat/_snaps/plots/priorpredictive.svg index 44d0414..bd09980 100644 --- a/tests/testthat/_snaps/plots/priorpredictive.svg +++ b/tests/testthat/_snaps/plots/priorpredictive.svg @@ -21,2272 +21,2267 @@ - - + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 @@ -2295,16 +2290,26 @@ - - - - -0 -200 -400 -600 -t + + + + + +0 +100 +200 +300 +400 +t me + +Censored + + + + +FALSE +TRUE priorpredictive diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R index 5746122..753d7f5 100644 --- a/tests/testthat/test-plots.R +++ b/tests/testthat/test-plots.R @@ -100,5 +100,25 @@ test_that("Can plot population trajectories with data", { mod <- biokinetics$new(data = data, covariate_formula = ~0 + infection_history) fit <- mod$fit() trajectories <- mod$simulate_population_trajectories() - vdiffr::expect_doppelganger("populationtrajectories_data", plot(trajectories, data = data)) + plot <- plot(trajectories, data = data) + expect_equal(length(plot$scales$scales), 1) + vdiffr::expect_doppelganger("populationtrajectories_data", plot) +}) + +test_that("Can plot population trajectories with log scale input data", { + local_mocked_bindings( + stan_package_model = mock_model, .package = "instantiate" + ) + # note that this is using a pre-fitted model with very few iterations, so the + # fits won't look very good + data <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) + data <- convert_log2_scale(data) + mod <- biokinetics$new(data = data, + covariate_formula = ~0 + infection_history, + scale = "log") + fit <- mod$fit() + trajectories <- mod$simulate_population_trajectories() + plot <- plot(trajectories, data = data) + expect_equal(length(plot$scales$scales), 0) + vdiffr::expect_doppelganger("populationtrajectories_logscale", plot) }) From 7b2daeccd41ace6641984a4fc9963fd2225fc305 Mon Sep 17 00:00:00 2001 From: "alex.hill@gmail.com" Date: Thu, 31 Oct 2024 16:20:19 +0000 Subject: [PATCH 12/12] factor our censoreed points logic into helper fn --- R/plot.R | 22 ++++++++++++---------- tests/testthat/test-plots.R | 2 +- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/R/plot.R b/R/plot.R index 6cc2a33..d8949ce 100644 --- a/R/plot.R +++ b/R/plot.R @@ -47,11 +47,7 @@ plot.biokinetics_priors <- function(x, geom_ribbon(aes(x = t, ymin = lo, ymax = hi), alpha = 0.5) if (!is.null(data)) { - if (!("censored" %in% colnames(data))) { - data$censored <- FALSE - } else { - data$censored <- data$censored != 0 - } + add_censored_indicator(data) dat <- data[time_since_last_exp <= tmax,] plot <- plot + geom_point(data = dat, size = 0.5, @@ -106,11 +102,7 @@ plot.biokinetics_population_trajectories <- function(x, ..., data = NULL) { } if (!is.null(data)) { validate_required_cols(data) - if (!("censored" %in% colnames(data))) { - data$censored <- FALSE - } else { - data$censored <- data$censored != 0 - } + add_censored_indicator(data) plot <- plot + geom_point(data = data, aes(x = as.integer(day - last_exp_day, units = "days"), @@ -129,3 +121,13 @@ plot.biokinetics_population_trajectories <- function(x, ..., data = NULL) { facet_formula <- function(covariates) { paste("~", paste(c("titre_type", covariates), collapse = "+")) } + +add_censored_indicator <- function(data) { + if (!("censored" %in% colnames(data))) { + # censored is an optional column in input data + # if not present, treat all points as uncensored + data[, censored:= FALSE] + } else { + data[, censored:= censored != 0] + } +} diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R index 753d7f5..710ec71 100644 --- a/tests/testthat/test-plots.R +++ b/tests/testthat/test-plots.R @@ -58,7 +58,7 @@ mock_model_multiple_covariates <- function(name, package) { list(sample = function(x, ...) readRDS(test_path("testdata", "testdraws_multiplecovariates.rds"))) } -test_that("Cna plot summarised and un-summarised population trajectories", { +test_that("Can plot summarised and un-summarised population trajectories", { # note that this is using a pre-fitted model with very few iterations, so the # fits won't look very good local_mocked_bindings(