condiitonal log transform and shape based on censoring

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

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

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

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

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