censor points based on limits

NAMESPACE

@@ -5,6 +5,7 @@ S3method(plot,biokinetics_priors)
 export(add_exposure_data)
 export(biokinetics)
 export(biokinetics_priors)
+export(convert_log2_scale)
 export(convert_log2_scale_inverse)
 export(plot_sero_data)
 importFrom(R6,R6Class)

R/cpp11.R

@@ -1,7 +1,7 @@
 # Generated by cpp11: do not edit by hand
 
-convert_log2_scale_inverse_cpp <- function(dt, vars_to_transform, lower_limit) {
-  .Call(`_epikinetics_convert_log2_scale_inverse_cpp`, dt, vars_to_transform, lower_limit)
+convert_log2_scale_inverse_cpp <- function(dt, vars_to_transform, smallest_value) {
+  .Call(`_epikinetics_convert_log2_scale_inverse_cpp`, dt, vars_to_transform, smallest_value)
 }
 
 simulate_trajectories_cpp <- function(person_params) {

R/plot.R

@@ -9,11 +9,15 @@
 #' @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.
+#' @param upper_detection_limit Optional upper detection limit.
+#' @param lower_detection_limit Optional lower detection limit.
 plot.biokinetics_priors <- function(x,
                                     ...,
                                     tmax = 150,
                                     n_draws = 2000,
-                                    data = NULL) {
+                                    data = NULL,
+                                    upper_detection_limit = NULL,
+                                    lower_detection_limit = NULL) {
 
   # Declare variables to suppress notes when compiling package
   # https://github.com/Rdatatable/data.table/issues/850#issuecomment-259466153
@@ -47,14 +51,13 @@ plot.biokinetics_priors <- function(x,
     geom_ribbon(aes(x = t, ymin = lo, ymax = hi), alpha = 0.5)
 
   if (!is.null(data)) {
-    add_censored_indicator(data)
     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"))
+                     y = value))
+
+    plot <- add_limits(plot, upper_detection_limit, lower_detection_limit)
   }
   plot
 }
@@ -67,51 +70,27 @@ plot.biokinetics_priors <- function(x,
 #' @param data A data.table with required columns time_since_last_exp, value and titre_type.
 #' @param tmax Integer. The number of time points in each simulated trajectory. Default 150.
 #' @param covariates Optional vector of covariate names to facet by (these must correspond to columns in the data.table)
-#' @param upper_detection_limit Optional upper detection limit. If provided, will be represented as a dashed line.
-#' @param lower_detection_limit Optional lower detection limit. If provided, will be represented as a dashed line.
+#' @param upper_detection_limit Optional upper detection limit.
+#' @param lower_detection_limit Optional lower detection limit.
 plot_sero_data <- function(data,
                            tmax = 150,
                            covariates = character(0),
                            upper_detection_limit = NULL,
                            lower_detection_limit = NULL) {
   validate_required_cols(data, c("time_since_last_exp", "value", "titre_type"))
-  add_censored_indicator(data)
-  data <- data[time_since_last_exp <= tmax, ]
+  data <- data[time_since_last_exp <= tmax,]
   # 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 <- censored <- NULL
 
   plot <- ggplot(data) +
-    geom_point(aes(x = time_since_last_exp, y = value, colour = titre_type, shape = censored),
+    geom_point(aes(x = time_since_last_exp, y = value, colour = titre_type),
                size = 0.5, alpha = 0.5) +
     geom_smooth(aes(x = time_since_last_exp, y = value, colour = titre_type)) +
     facet_wrap(eval(parse(text = facet_formula(covariates)))) +
-    guides(shape = guide_legend(title = "Censored"),
-           colour = guide_legend(title = "Titre type"))
+    guides(colour = guide_legend(title = "Titre type"))
 
-  if (!is.null(lower_detection_limit)) {
-    plot <- plot +
-      geom_hline(yintercept = lower_detection_limit,
-                 linetype = 'dotted') +
-      annotate("text", x = 1,
-               y = lower_detection_limit,
-               label = "Lower detection limit",
-               vjust = -0.5,
-               hjust = 0,
-               size = 3)
-  }
-  if (!is.null(upper_detection_limit)) {
-    plot <- plot +
-      geom_hline(yintercept = upper_detection_limit,
-                 linetype = 'dotted') +
-      annotate("text", x = 1,
-               y = upper_detection_limit,
-               label = "Upper detection limit",
-               vjust = -0.5,
-               hjust = 0,
-               size = 3)
-  }
-  plot
+  add_limits(plot, upper_detection_limit, lower_detection_limit)
 }
 
 #' Plot method for "biokinetics_population_trajectories" class
@@ -122,13 +101,16 @@ plot_sero_data <- function(data,
 #' @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, ..., data = NULL) {
+plot.biokinetics_population_trajectories <- function(x, ...,
+                                                     data = NULL) {
   covariates <- attr(x, "covariates")
+  upper_detection_limit <- attr(x, "upper_detection_limit")
+  lower_detection_limit <- attr(x, "lower_detection_limit")
 
   # 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 <- censored <- mu <- .draw <- NULL
+  day <- last_exp_day <- mu <- .draw <- NULL
 
   if (attr(x, "summarised")) {
     plot <- ggplot(x) +
@@ -144,35 +126,50 @@ plot.biokinetics_population_trajectories <- function(x, ..., data = NULL) {
   }
   if (!is.null(data)) {
     validate_required_cols(data)
-    add_censored_indicator(data)
     plot <- plot +
       geom_point(data = data,
                  aes(x = as.integer(day - last_exp_day, units = "days"),
-                     y = value, shape = censored), size = 0.5, alpha = 0.5)
+                     y = value), size = 0.5, alpha = 0.5)
   }
   if (attr(x, "scale") == "natural") {
     plot <- plot + scale_y_continuous(trans = "log2")
   }
-  plot +
+  plot <- plot +
     facet_wrap(eval(parse(text = facet_formula(covariates)))) +
     guides(fill = guide_legend(title = "Titre type"),
-           colour = "none",
-           shape = guide_legend(title = "Censored"))
+           colour = "none")
+  if (!is.null(data)) {
+    plot <- add_limits(plot, upper_detection_limit, lower_detection_limit)
+  }
+  plot
 }
 
 facet_formula <- function(covariates) {
   paste("~", paste(c("titre_type", covariates), collapse = "+"))
 }
 
-add_censored_indicator <- function(data) {
-  # Declare variables to suppress notes when compiling package
-  # https://github.com/Rdatatable/data.table/issues/850#issuecomment-259466153
-  censored <- NULL
-  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]
+add_limits <- function(plot, upper_detection_limit, lower_detection_limit) {
+  if (!is.null(lower_detection_limit)) {
+    plot <- plot +
+      geom_hline(yintercept = lower_detection_limit,
+                 linetype = 'dotted') +
+      annotate("text", x = 1,
+               y = lower_detection_limit,
+               label = "Lower detection limit",
+               vjust = -0.5,
+               hjust = 0,
+               size = 3)
   }
-}
+  if (!is.null(upper_detection_limit)) {
+    plot <- plot +
+      geom_hline(yintercept = upper_detection_limit,
+                 linetype = 'dotted') +
+      annotate("text", x = 1,
+               y = upper_detection_limit,
+               label = "Upper detection limit",
+               vjust = -0.5,
+               hjust = 0,
+               size = 3)
+  }
+  plot
+}

R/utils.R

@@ -1,10 +1,20 @@
+#' @title Base 2 log scale conversion
+#'
+#' @description Natural scale data is converted to a base 2 log scale before model fitting.
+#' This function does not modify the provided data.table in-place, but returns a transformed copy.
+#' @return A data.table, identical to the input data but with specified columns transformed.
+#' @param dt_in data.table containing data to be transformed to base 2 log scale.
+#' @param vars_to_transform Names of columns to apply the transformation to.
+#' @param smallest_value The lowest value in the original dataset. Used to get the data onto
+#' a scale that starts at zero.
+#' @export
 convert_log2_scale <- function(
   dt_in,
-  lower_limit,
+  smallest_value,
   vars_to_transform = "value") {
   dt_out <- data.table::copy(dt_in)
   for (var in vars_to_transform) {
-    dt_out[, (var) := log2(get(var) / lower_limit)]
+    dt_out[, (var) := log2(get(var) / smallest_value)]
   }
   return(dt_out)
 }
@@ -17,12 +27,12 @@ convert_log2_scale <- function(
 #' @return A data.table, identical to the input data but with specified columns transformed.
 #' @param dt_in data.table containing data to be transformed from base 2 log to natural scale.
 #' @param vars_to_transform Names of columns to apply the transformation to.
-#' @param lower_limit The lowest value in the original dataset.
+#' @param smallest_value The lowest value in the original dataset.
 #' @export
-convert_log2_scale_inverse <- function(dt_in, vars_to_transform, lower_limit) {
+convert_log2_scale_inverse <- function(dt_in, vars_to_transform, smallest_value) {
   dt_out <- data.table::copy(dt_in)
   for (var in vars_to_transform) {
-    dt_out[, (var) := lower_limit * 2^(get(var))]
+    dt_out[, (var) := smallest_value * 2^(get(var))]
   }
   dt_out
 }