trajectory plots and doppelganger tests

DESCRIPTION

@@ -30,7 +30,7 @@ Suggests:
     knitr,
     lubridate,
     rmarkdown,
-    testthat
+    testthat,
     vdiffr
 VignetteBuilder: knitr
 LinkingTo: 

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)

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
     }
   )
 )

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
 

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 = "+"))
+}

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)