From 730a4218557af4bbeefb3eadf4be506aae41edfd Mon Sep 17 00:00:00 2001 From: "alex.hill@gmail.com" Date: Thu, 14 Nov 2024 15:44:52 +0000 Subject: [PATCH] add individual trajectories plot --- R/biokinetics.R | 2 +- R/plot.R | 59 +++++++++++- .../_snaps/plots/individualtrajectories.svg | 96 +++++++++++++++++++ tests/testthat/test-plots.R | 16 ++++ 4 files changed, 170 insertions(+), 3 deletions(-) create mode 100644 tests/testthat/_snaps/plots/individualtrajectories.svg diff --git a/R/biokinetics.R b/R/biokinetics.R index d7f57eb..a328d0a 100644 --- a/R/biokinetics.R +++ b/R/biokinetics.R @@ -475,7 +475,7 @@ biokinetics <- R6::R6Class( }, #' @description Simulate individual trajectories from the model. This is #' computationally expensive and may take a while to run if n_draws is large. - #' @return A data.table. If summarise = TRUE columns are calendar_date, titre_type, me, lo, hi, time_shift. + #' @return A data.table. If summarise = TRUE columns are calendar_day, titre_type, me, lo, hi, time_shift. #' If summarise = FALSE, columns are pid, draw, time_since_last_exp, mu, titre_type, exposure_day, calendar_day, time_shift #' and a column for each covariate in the regression model. See the data vignette for details: #' \code{vignette("data", package = "epikinetics")}. diff --git a/R/plot.R b/R/plot.R index d8949ce..a00668d 100644 --- a/R/plot.R +++ b/R/plot.R @@ -118,6 +118,61 @@ plot.biokinetics_population_trajectories <- function(x, ..., data = NULL) { shape = guide_legend(title = "Censored")) } +#' Plot method for "biokinetics_individual_trajectories" class +#' +#' @param x An object of class "biokinetics_individual_trajectories". These are +#' generated by running biokinetics$simulate_individual_trajectories(). See +#' \href{../../epikinetics/html/biokinetics.html#method-biokinetics-simulate_individaul_trajectories}{\code{biokinetics$simulate_individual_trajectories()}} +#' @param \dots Further arguments passed to the method. +#' @export +plot.biokinetics_individual_trajectories <- function(x, ..., data = NULL, + min_date = NULL, + max_date = NULL) { + + # Declare variables to suppress notes when compiling package + # https://github.com/Rdatatable/data.table/issues/850#issuecomment-259466153 + calendar_day <- value <- me <- mu <- titre_type <- lo <- hi <- day <- pid <- NULL + if (is.null(min_date)) { + min_date <- min(x$calendar_day) + } + if (is.null(max_date)) { + max_date <- max(x$calendar_day) + } + if (attr(x, "summarised")) { + plot <- ggplot(x) + + geom_line(aes(x = calendar_day, y = me, group = titre_type, colour = titre_type)) + + geom_ribbon(aes(x = calendar_day, + ymin = lo, + ymax = hi, + group = titre_type), alpha = 0.5) + } else { + plot <- ggplot(x) + + geom_line(aes(x = calendar_day, y = mu, + colour = titre_type, group = pid), alpha = 0.1, linewidth = 0.1) + } + if (!is.null(data)) { + validate_required_cols(data, c("day", "value")) + plot <- plot + + geom_point(data = data, + aes(x = day, + y = value), size = 0.5, alpha = 0.5) + } + plot + + labs(x = "Date", + y = expression(paste("Titre (IC"[50], ")"))) + + geom_smooth( + aes(x = calendar_day, + y = me, + fill = titre_type, + colour = titre_type, + group = titre_type), + alpha = 0.5, span = 0.2) + + scale_x_date(date_labels = "%b %Y", + limits = c(min_date, max_date)) + + guides(colour = guide_legend(title = "Titre type"), + fill = "none") +} + facet_formula <- function(covariates) { paste("~", paste(c("titre_type", covariates), collapse = "+")) } @@ -126,8 +181,8 @@ 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] + data[, censored := FALSE] } else { - data[, censored:= censored != 0] + data[, censored := censored != 0] } } diff --git a/tests/testthat/_snaps/plots/individualtrajectories.svg b/tests/testthat/_snaps/plots/individualtrajectories.svg new file mode 100644 index 0000000..2b3a4c2 --- /dev/null +++ b/tests/testthat/_snaps/plots/individualtrajectories.svg @@ -0,0 +1,96 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +0 +300 +600 +900 + + + + + + + + +Jan 2021 +Jul 2021 +Jan 2022 +Jul 2022 +Date +Titre (IC +50 +) + +Titre type + + + + + + + + + + + + +Alpha +Ancestral +Delta +individualtrajectories + + diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R index ccb2bcf..6131e52 100644 --- a/tests/testthat/test-plots.R +++ b/tests/testthat/test-plots.R @@ -122,3 +122,19 @@ test_that("Can plot population trajectories with log scale input data", { expect_equal(length(plot$scales$scales), 0) vdiffr::expect_doppelganger("populationtrajectories_logscale", plot) }) + +test_that("Can plot summarised individual 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( + stan_package_model = mock_model, .package = "instantiate" + ) + mod <- biokinetics$new(file_path = system.file("delta_full.rds", package = "epikinetics"),) + mod$fit() + trajectories <- mod$simulate_individual_trajectories(n_draws = 250, + summarise = TRUE) + # because these fits are so bad there are some v high upper values, so just + # create these articially + trajectories[, hi := me + 100] + vdiffr::expect_doppelganger("individualtrajectories", plot(trajectories)) +})