diff --git a/R/biokinetics.R b/R/biokinetics.R index c6711b2..d087f02 100644 --- a/R/biokinetics.R +++ b/R/biokinetics.R @@ -11,8 +11,8 @@ biokinetics <- R6::R6Class( scale = NULL, priors = NULL, preds_sd = NULL, - upper_detection_limit = NULL, - lower_detection_limit = NULL, + upper_censoring_limit = NULL, + lower_censoring_limit = NULL, smallest_value = NULL, data = NULL, covariate_formula = NULL, @@ -29,28 +29,18 @@ biokinetics <- R6::R6Class( stop("Model has not been fitted yet. Call 'fit' before calling this function.") } }, - get_upper_detection_limit = function(upper_detection_limit) { - if (is.null(upper_detection_limit)) { - private$upper_detection_limit <- max(private$data$value) + get_upper_censoring_limit = function(upper_censoring_limit) { + if (is.null(upper_censoring_limit)) { + private$upper_censoring_limit <- max(private$data$value) } else { - max_value <- max(private$data$value) - if (max_value >= upper_detection_limit) { - warning(sprintf("Data contains values >= the upper detection limit %s. These will be censored.", - upper_detection_limit)) - } - private$upper_detection_limit <- upper_detection_limit + private$upper_censoring_limit <- upper_censoring_limit } }, - get_lower_detection_limit = function(lower_detection_limit) { - if (is.null(lower_detection_limit)) { - private$lower_detection_limit <- min(private$data$value) + get_lower_censoring_limit = function(lower_censoring_limit) { + if (is.null(lower_censoring_limit)) { + private$lower_censoring_limit <- min(private$data$value) } else { - min_value <- min(private$data$value) - if (min_value <= lower_detection_limit) { - warning(sprintf("Data contains values <= the lower detection limit %s. These will be censored.", - lower_detection_limit)) - } - private$lower_detection_limit <- lower_detection_limit + private$lower_censoring_limit <- lower_censoring_limit } }, model_matrix_with_dummy = function(data) { @@ -204,11 +194,11 @@ biokinetics <- R6::R6Class( stan_data$P <- ncol(private$design_matrix) if (private$scale == "natural") { # do the same transformation as used on the data - stan_data$upper_detection_limit <- log2(private$upper_detection_limit / private$smallest_value) - stan_data$lower_detection_limit <- log2(private$lower_detection_limit / private$smallest_value) + stan_data$upper_censoring_limit <- log2(private$upper_censoring_limit / private$smallest_value) + stan_data$lower_censoring_limit <- log2(private$lower_censoring_limit / private$smallest_value) } else { - stan_data$upper_detection_limit <- private$upper_detection_limit - stan_data$lower_detection_limit <- private$lower_detection_limit + stan_data$upper_censoring_limit <- private$upper_censoring_limit + stan_data$lower_censoring_limit <- private$lower_censoring_limit } private$stan_input_data <- c(stan_data, private$priors) }, @@ -248,24 +238,24 @@ biokinetics <- R6::R6Class( #' @param preds_sd Standard deviation of predictor coefficients. Default 0.25. #' @param scale One of "log" or "natural". Default "natural". Is provided data on a log or a natural scale? If on a natural scale it #' will be converted to a log scale for model fitting. - #' @param upper_detection_limit Optional upper detection limit of the titre used. This is needed to construct a likelihood for upper censored - #' values, so only needs to be provided if you have such values in the dataset. If not provided, the model will default - #' to using the largest value in the dataset as the upper detection limit. - #' @param lower_detection_limit Optional lower detection limit of the titre used. This is needed to construct a likelihood for lower censored - #' values, so only needs to be provided if you have such values in the dataset. If not provided, the model will default - #' to using the smallest value in the dataset as the lower detection limit. - #' @param truncate_upper Logical. Whether values greater than the upper detection limit should be truncated, i.e. set to the same value as the limit. Default TRUE. - #' @param truncate_lower Logical. Whether values smaller than the lower detection limit should be truncated, i.e. set to the same value as the limit. Default TRUE. + #' @param upper_censoring_limit Optional value at which to upper censor data points. This is needed to construct a likelihood for upper censored + #' values, so only needs to be provided if you have such values in the dataset. If not provided, no censoring will be done. + #' @param lower_censoring_limit Optional value at which to lower censor data points. This is needed to construct a likelihood for lower censored + #' values, so only needs to be provided if you have such values in the dataset. If not provided, no censoring will be done. + #' @param strict_upper_limit Logical. Whether values greater than the upper censoring limit should be censored. + #' If FALSE, only values exactly equal to the upper censoring limit will be censored. Default TRUE. + #' @param strict_lower_limit Logical. Whether values smaller than the lower censoring limit should be censored. + #' If FALSE, only values exactly equal to the lower censoring limit will be censored. Default TRUE. initialize = function(priors = biokinetics_priors(), data = NULL, file_path = NULL, covariate_formula = ~0, preds_sd = 0.25, scale = "natural", - upper_detection_limit = NULL, - lower_detection_limit = NULL, - truncate_upper = TRUE, - truncate_lower = TRUE) { + upper_censoring_limit = NULL, + lower_censoring_limit = NULL, + strict_upper_limit = TRUE, + strict_lower_limit = TRUE) { validate_priors(priors) private$priors <- priors validate_numeric(preds_sd) @@ -291,18 +281,36 @@ biokinetics <- R6::R6Class( validate_required_cols(private$data) validate_formula_vars(private$all_formula_vars, private$data) logger::log_info("Preparing data for stan") - private$get_upper_detection_limit(upper_detection_limit) - private$get_lower_detection_limit(lower_detection_limit) - if (truncate_upper) { - private$data[, value := ifelse(value > private$upper_detection_limit, private$upper_detection_limit, value)] - } - if (truncate_lower) { - private$data[, value := ifelse(value < private$lower_detection_limit, private$lower_detection_limit, value)] + private$get_upper_censoring_limit(upper_censoring_limit) + private$get_lower_censoring_limit(lower_censoring_limit) + max_value <- max(private$data$value) + min_value <- min(private$data$value) + values_above <- max_value > private$upper_censoring_limit + values_below <- min_value < private$lower_censoring_limit + if (strict_upper_limit) { + if (values_above) { + warning(sprintf("Data contains values above the upper censoring limit %s and these will be censored. To turn off this behaviour set strict_upper_limit to FALSE.", + private$upper_censoring_limit)) + } + private$data[, value := ifelse(value > private$upper_censoring_limit, private$upper_censoring_limit, value)] + } else if (values_above) { + warning(sprintf("Data contains values above the upper censoring limit %s. To treat these as censored set strict_upper_limit to TRUE.", + private$upper_censoring_limit)) + } + if (strict_lower_limit) { + if (values_below) { + warning(sprintf("Data contains values below the lower censoring limit %s and these will be censored. To turn off this behaviour set strict_lower_limit to FALSE.", + private$lower_censoring_limit)) + } + private$data[, value := ifelse(value < private$lower_censoring_limit, private$lower_censoring_limit, value)] + } else if (values_below) { + warning(sprintf("Data contains values below the lower censoring limit %s. To treat these as censored set strict_lower_limit to TRUE.", + private$lower_censoring_limit)) } private$data[, `:=`(obs_id = seq_len(.N), time_since_last_exp = as.integer(day - last_exp_day, units = "days"), - censored_lo = value == private$lower_detection_limit, - censored_hi = value == private$upper_detection_limit)] + censored_lo = value == private$lower_censoring_limit, + censored_hi = value == private$upper_censoring_limit)] private$construct_design_matrix() private$build_covariate_lookup_table() private$build_pid_lookup() @@ -332,8 +340,8 @@ biokinetics <- R6::R6Class( tmax = tmax, n_draws = n_draws, data = private$data, - upper_detection_limit = private$stan_input_data$upper_detection_limit, - lower_detection_limit = private$stan_input_data$lower_detection_limit) + upper_censoring_limit = private$stan_input_data$upper_censoring_limit, + lower_censoring_limit = private$stan_input_data$lower_censoring_limit) }, #' @description Plot model input data with a smoothing function. Note that #' this plot is of the data as provided to the Stan model so is on a log scale, @@ -344,32 +352,32 @@ biokinetics <- R6::R6Class( plot_sero_data(private$data, tmax = tmax, covariates = private$all_formula_vars, - upper_detection_limit = private$stan_input_data$upper_detection_limit, - lower_detection_limit = private$stan_input_data$lower_detection_limit) + upper_censoring_limit = private$stan_input_data$upper_censoring_limit, + lower_censoring_limit = private$stan_input_data$lower_censoring_limit) }, - #' @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. + #' @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() { private$stan_input_data }, - #' @description View the mapping of human readable covariate names to the model variable p. - #' @return A data.table mapping the model variable p to human readable covariates. + #' @description View the mapping of human readable covariate names to the model variable p. + #' @return A data.table mapping the model variable p to human readable covariates. get_covariate_lookup_table = function() { private$covariate_lookup_table }, - #' @description Fit the model and return CmdStanMCMC fitted model object. - #' @return A CmdStanMCMC fitted model object: - #' @param ... Named arguments to the `sample()` method of CmdStan model. - #' objects: + #' @description Fit the model and return CmdStanMCMC fitted model object. + #' @return A CmdStanMCMC fitted model object: + #' @param ... Named arguments to the `sample()` method of CmdStan model. + #' objects: fit = function(...) { logger::log_info("Fitting model") private$fitted <- private$model$sample(private$stan_input_data, ...) private$fitted }, - #' @description Extract fitted population parameters - #' @return A data.table - #' @param n_draws Integer. Default 2000. - #' @param human_readable_covariates Logical. Default TRUE. + #' @description Extract fitted population parameters + #' @return A data.table + #' @param n_draws Integer. Default 2000. + #' @param human_readable_covariates Logical. Default TRUE. extract_population_parameters = function(n_draws = 2000, human_readable_covariates = TRUE) { private$check_fitted() @@ -489,8 +497,8 @@ biokinetics <- R6::R6Class( attr(dt_out, "summarised") <- summarise attr(dt_out, "scale") <- private$scale attr(dt_out, "covariates") <- private$all_formula_vars - attr(dt_out, "upper_detection_limit") <- private$upper_detection_limit - attr(dt_out, "lower_detection_limit") <- private$lower_detection_limit + attr(dt_out, "upper_censoring_limit") <- private$upper_censoring_limit + attr(dt_out, "lower_censoring_limit") <- private$lower_censoring_limit dt_out }, diff --git a/R/plot.R b/R/plot.R index 56346eb..cb49ca8 100644 --- a/R/plot.R +++ b/R/plot.R @@ -9,15 +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. +#' @param upper_censoring_limit Optional upper detection limit. +#' @param lower_censoring_limit Optional lower detection limit. plot.biokinetics_priors <- function(x, ..., tmax = 150, n_draws = 2000, data = NULL, - upper_detection_limit = NULL, - lower_detection_limit = NULL) { + upper_censoring_limit = NULL, + lower_censoring_limit = NULL) { # Declare variables to suppress notes when compiling package # https://github.com/Rdatatable/data.table/issues/850#issuecomment-259466153 @@ -57,7 +57,7 @@ plot.biokinetics_priors <- function(x, aes(x = time_since_last_exp, y = value)) - plot <- add_limits(plot, upper_detection_limit, lower_detection_limit) + plot <- add_limits(plot, upper_censoring_limit, lower_censoring_limit) } plot } @@ -70,13 +70,13 @@ 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. -#' @param lower_detection_limit Optional lower detection limit. +#' @param upper_censoring_limit Optional upper detection limit. +#' @param lower_censoring_limit Optional lower detection limit. plot_sero_data <- function(data, tmax = 150, covariates = character(0), - upper_detection_limit = NULL, - lower_detection_limit = NULL) { + upper_censoring_limit = NULL, + lower_censoring_limit = NULL) { validate_required_cols(data, c("time_since_last_exp", "value", "titre_type")) data <- data[time_since_last_exp <= tmax,] # Declare variables to suppress notes when compiling package @@ -90,7 +90,7 @@ plot_sero_data <- function(data, facet_wrap(eval(parse(text = facet_formula(covariates)))) + guides(colour = guide_legend(title = "Titre type")) - add_limits(plot, upper_detection_limit, lower_detection_limit) + add_limits(plot, upper_censoring_limit, lower_censoring_limit) } #' Plot method for "biokinetics_population_trajectories" class @@ -104,8 +104,8 @@ plot_sero_data <- function(data, 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") + upper_censoring_limit <- attr(x, "upper_censoring_limit") + lower_censoring_limit <- attr(x, "lower_censoring_limit") # Declare variables to suppress notes when compiling package # https://github.com/Rdatatable/data.table/issues/850#issuecomment-259466153 @@ -139,7 +139,7 @@ plot.biokinetics_population_trajectories <- function(x, ..., guides(fill = guide_legend(title = "Titre type"), colour = "none") if (!is.null(data)) { - plot <- add_limits(plot, upper_detection_limit, lower_detection_limit) + plot <- add_limits(plot, upper_censoring_limit, lower_censoring_limit) } plot } @@ -148,24 +148,24 @@ facet_formula <- function(covariates) { paste("~", paste(c("titre_type", covariates), collapse = "+")) } -add_limits <- function(plot, upper_detection_limit, lower_detection_limit) { - if (!is.null(lower_detection_limit)) { +add_limits <- function(plot, upper_censoring_limit, lower_censoring_limit) { + if (!is.null(lower_censoring_limit)) { plot <- plot + - geom_hline(yintercept = lower_detection_limit, + geom_hline(yintercept = lower_censoring_limit, linetype = 'dotted') + annotate("text", x = 1, - y = lower_detection_limit, + y = lower_censoring_limit, label = "Lower detection limit", vjust = -0.5, hjust = 0, size = 3) } - if (!is.null(upper_detection_limit)) { + if (!is.null(upper_censoring_limit)) { plot <- plot + - geom_hline(yintercept = upper_detection_limit, + geom_hline(yintercept = upper_censoring_limit, linetype = 'dotted') + annotate("text", x = 1, - y = upper_detection_limit, + y = upper_censoring_limit, label = "Upper detection limit", vjust = -0.5, hjust = 0, diff --git a/man/biokinetics.Rd b/man/biokinetics.Rd index 7e3b79f..7d6a65c 100644 --- a/man/biokinetics.Rd +++ b/man/biokinetics.Rd @@ -36,10 +36,10 @@ Initialise the kinetics model. covariate_formula = ~0, preds_sd = 0.25, scale = "natural", - upper_detection_limit = NULL, - lower_detection_limit = NULL, - truncate_upper = TRUE, - truncate_lower = TRUE + upper_censoring_limit = NULL, + lower_censoring_limit = NULL, + strict_upper_limit = TRUE, + strict_lower_limit = TRUE )}\if{html}{\out{}} } @@ -61,17 +61,17 @@ will be treated as categorical variables. Default ~0.} \item{\code{scale}}{One of "log" or "natural". Default "natural". Is provided data on a log or a natural scale? If on a natural scale it will be converted to a log scale for model fitting.} -\item{\code{upper_detection_limit}}{Optional upper detection limit of the titre used. This is needed to construct a likelihood for upper censored -values, so only needs to be provided if you have such values in the dataset. If not provided, the model will default -to using the largest value in the dataset as the upper detection limit.} +\item{\code{upper_censoring_limit}}{Optional value at which to upper censor data points. This is needed to construct a likelihood for upper censored +values, so only needs to be provided if you have such values in the dataset. If not provided, no censoring will be done.} -\item{\code{lower_detection_limit}}{Optional lower detection limit of the titre used. This is needed to construct a likelihood for lower censored -values, so only needs to be provided if you have such values in the dataset. If not provided, the model will default -to using the smallest value in the dataset as the lower detection limit.} +\item{\code{lower_censoring_limit}}{Optional value at which to lower censor data points. This is needed to construct a likelihood for lower censored +values, so only needs to be provided if you have such values in the dataset. If not provided, no censoring will be done.} -\item{\code{truncate_upper}}{Logical. Whether values greater than the upper detection limit should be truncated, i.e. set to the same value as the limit. Default TRUE.} +\item{\code{strict_upper_limit}}{Logical. Whether values greater than the upper censoring limit should be censored. +If FALSE, only values exactly equal to the upper censoring limit will be censored. Default TRUE.} -\item{\code{truncate_lower}}{Logical. Whether values smaller than the lower detection limit should be truncated, i.e. set to the same value as the limit. Default TRUE.} +\item{\code{strict_lower_limit}}{Logical. Whether values smaller than the lower censoring limit should be censored. +If FALSE, only values exactly equal to the lower censoring limit will be censored. Default TRUE.} } \if{html}{\out{}} } diff --git a/man/plot.biokinetics_priors.Rd b/man/plot.biokinetics_priors.Rd index 56cb964..4db8eba 100644 --- a/man/plot.biokinetics_priors.Rd +++ b/man/plot.biokinetics_priors.Rd @@ -11,8 +11,8 @@ and optionally compare to a dataset.} tmax = 150, n_draws = 2000, data = NULL, - upper_detection_limit = NULL, - lower_detection_limit = NULL + upper_censoring_limit = NULL, + lower_censoring_limit = NULL ) } \arguments{ @@ -26,9 +26,9 @@ and optionally compare to a dataset.} \item{data}{Optional data.frame with columns time_since_last_exp and value. The raw data to compare to.} -\item{upper_detection_limit}{Optional upper detection limit.} +\item{upper_censoring_limit}{Optional upper detection limit.} -\item{lower_detection_limit}{Optional lower detection limit.} +\item{lower_censoring_limit}{Optional lower detection limit.} } \value{ A ggplot2 object. diff --git a/man/plot_sero_data.Rd b/man/plot_sero_data.Rd index ec83395..6986bc0 100644 --- a/man/plot_sero_data.Rd +++ b/man/plot_sero_data.Rd @@ -8,8 +8,8 @@ plot_sero_data( data, tmax = 150, covariates = character(0), - upper_detection_limit = NULL, - lower_detection_limit = NULL + upper_censoring_limit = NULL, + lower_censoring_limit = NULL ) } \arguments{ @@ -19,9 +19,9 @@ plot_sero_data( \item{covariates}{Optional vector of covariate names to facet by (these must correspond to columns in the data.table)} -\item{upper_detection_limit}{Optional upper detection limit.} +\item{upper_censoring_limit}{Optional upper detection limit.} -\item{lower_detection_limit}{Optional lower detection limit.} +\item{lower_censoring_limit}{Optional lower detection limit.} } \value{ A ggplot2 object. diff --git a/src/stan/antibody_kinetics_main.stan b/src/stan/antibody_kinetics_main.stan index 48b1e4e..2d6cb38 100644 --- a/src/stan/antibody_kinetics_main.stan +++ b/src/stan/antibody_kinetics_main.stan @@ -29,8 +29,8 @@ data { int N; // Number of observations int N_events; // Number of exposure events int K; // Number of titre types - real upper_detection_limit; // Upper detection limit - real lower_detection_limit; // Lower detection limit + real upper_censoring_limit; // Upper detection limit + real lower_censoring_limit; // Lower detection limit array[N] int titre_type; // Titre type for each observation vector[N] t; // Time for each observation vector[N] value; // Observed titre values @@ -148,10 +148,10 @@ model { value[uncens_idx] ~ normal(mu[uncens_idx], sigma); // Likelihood for observations at lower limit - target += normal_lcdf(lower_detection_limit | mu[cens_lo_idx], sigma); + target += normal_lcdf(lower_censoring_limit | mu[cens_lo_idx], sigma); // Censoring at upper limit - target += normal_lccdf(upper_detection_limit | mu[cens_hi_idx], sigma); + target += normal_lccdf(upper_censoring_limit | mu[cens_hi_idx], sigma); // Covariate-level mean priors, parameterised from previous studies t0_pop ~ normal(mu_t0, sigma_t0); diff --git a/tests/testthat/test-data.R b/tests/testthat/test-data.R index 2a7f005..4878b00 100644 --- a/tests/testthat/test-data.R +++ b/tests/testthat/test-data.R @@ -19,7 +19,7 @@ test_that("Can construct stan data", { expect_equal(names(stan_data), c("N", "N_events", "id", "value", "titre_type", "preds_sd", "K", "N_uncens", "N_lo", "N_hi", "uncens_idx", "cens_lo_idx", - "cens_hi_idx", "t", "X", "P", "upper_detection_limit", "lower_detection_limit", "mu_t0", + "cens_hi_idx", "t", "X", "P", "upper_censoring_limit", "lower_censoring_limit", "mu_t0", "mu_tp", "mu_ts", "mu_m1", "mu_m2", "mu_m3", "sigma_t0", "sigma_tp", "sigma_ts", "sigma_m1", "sigma_m2", "sigma_m3")) @@ -32,8 +32,8 @@ test_that("Can construct stan data", { test_that("Data above/below limits is censored", { dat <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) suppressWarnings({ mod <- biokinetics$new(data = dat, - lower_detection_limit = 10, - upper_detection_limit = 2500) }) + lower_censoring_limit = 10, + upper_censoring_limit = 2500) }) stan_data <- mod$get_stan_data() expect_equal(stan_data$N_uncens, dat[value > 10 & value < 2500, .N]) expect_equal(stan_data$N_lo, dat[value <= 10, .N]) @@ -72,47 +72,68 @@ test_that("Log scale data is passed directly to stan", { test_that("Highest value is used as default upper limit", { dat <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) - mod <- biokinetics$new(data = dat, lower_detection_limit = 2) + mod <- biokinetics$new(data = dat, lower_censoring_limit = 2) stan_data <- mod$get_stan_data() - expect_equal(stan_data$upper_detection_limit, log2(max(dat$value) / min(dat$value))) - expect_equal(stan_data$lower_detection_limit, log2(2 / min(dat$value))) + expect_equal(stan_data$upper_censoring_limit, log2(max(dat$value) / min(dat$value))) + expect_equal(stan_data$lower_censoring_limit, log2(2 / min(dat$value))) }) -test_that("Warns if data contains values above the upper limit", { +test_that("Warns if data contains values above the upper limit and limit is strict", { dat <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) - expect_warning({ mod <- biokinetics$new(data = dat, upper_detection_limit = 10) }, - "Data contains values >= the upper detection limit 10. These will be censored.") + expect_warning({ mod <- biokinetics$new(data = dat, upper_censoring_limit = 10) }, + "Data contains values above the upper censoring limit 10 and these will be censored. To turn off this behaviour set strict_upper_limit to FALSE.") stan_data <- mod$get_stan_data() - expect_equal(stan_data$upper_detection_limit, log2(10 / min(dat$value))) - expect_equal(stan_data$lower_detection_limit, 0) + expect_equal(stan_data$upper_censoring_limit, log2(10 / min(dat$value))) + expect_equal(stan_data$lower_censoring_limit, 0) + expect_true(all(stan_data$value <= stan_data$upper_censoring_limit)) +}) + +test_that("Warns if data contains values above the upper limit and limit is not strict", { + dat <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) + expect_warning({ mod <- biokinetics$new(data = dat, upper_censoring_limit = 10, strict_upper_limit = FALSE) }, + "Data contains values above the upper censoring limit 10. To treat these as censored set strict_upper_limit to TRUE.") + stan_data <- mod$get_stan_data() + expect_equal(stan_data$upper_censoring_limit, log2(10 / min(dat$value))) + expect_equal(stan_data$lower_censoring_limit, 0) + expect_false(all(stan_data$value <= stan_data$upper_censoring_limit)) }) test_that("Smallest value is used as default lower limit", { dat <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) - mod <- biokinetics$new(data = dat, upper_detection_limit = 3000) + mod <- biokinetics$new(data = dat, upper_censoring_limit = 3000) stan_data <- mod$get_stan_data() - expect_equal(stan_data$upper_detection_limit, log2(3000 / min(dat$value))) - expect_equal(stan_data$lower_detection_limit, 0) + expect_equal(stan_data$upper_censoring_limit, log2(3000 / min(dat$value))) + expect_equal(stan_data$lower_censoring_limit, 0) }) test_that("Warns if data contains values below the lower limit", { dat <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) expect_warning({ mod <- biokinetics$new(data = dat, - lower_detection_limit = 10) }, - "Data contains values <= the lower detection limit 10. These will be censored.") + lower_censoring_limit = 10) }, + "Data contains values below the lower censoring limit 10 and these will be censored. To turn off this behaviour set strict_lower_limit to FALSE.") + stan_data <- mod$get_stan_data() + expect_equal(stan_data$lower_censoring_limit, log2(10/min(dat$value))) + expect_true(all(stan_data$value >= stan_data$lower_censoring_limit)) +}) + +test_that("Warns if data contains values below the lower limit and limit is not strict", { + dat <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) + expect_warning({ mod <- biokinetics$new(data = dat, lower_censoring_limit = 10, strict_lower_limit = FALSE) }, + "Data contains values below the lower censoring limit 10. To treat these as censored set strict_lower_limit to TRUE.") stan_data <- mod$get_stan_data() - expect_equal(stan_data$lower_detection_limit, log2(10/min(dat$value))) + expect_equal(stan_data$lower_censoring_limit, log2(10/min(dat$value))) + expect_false(all(stan_data$value <= stan_data$lower_censoring_limit)) }) -test_that("Detection limits are passed to stan without transformation if using log data", { +test_that("Censoring limits are passed to stan without transformation if using log data", { dat <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) mod <- biokinetics$new(data = convert_log2_scale(dat, smallest_value = 2, vars_to_transform = "value"), scale = "log", - lower_detection_limit = 1, - upper_detection_limit = 15) + lower_censoring_limit = 1, + upper_censoring_limit = 15) stan_data <- mod$get_stan_data() - expect_equal(stan_data$upper_detection_limit, 15) - expect_equal(stan_data$lower_detection_limit, 1) + expect_equal(stan_data$upper_censoring_limit, 15) + expect_equal(stan_data$lower_censoring_limit, 1) }) diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R index ec00126..1808e51 100644 --- a/tests/testthat/test-plots.R +++ b/tests/testthat/test-plots.R @@ -42,15 +42,15 @@ test_that("Prior predictions from model are the same", { test_that("Can plot input data", { data <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) mod <- biokinetics$new(data = data, - lower_detection_limit = 40, - truncate_lower = FALSE) + lower_censoring_limit = 40, + strict_lower_limit = FALSE) suppressWarnings({plot <- mod$plot_model_inputs()}) vdiffr::expect_doppelganger("inputdata", plot) mod <- biokinetics$new(data = data, - lower_detection_limit = 40, + lower_censoring_limit = 40, covariate_formula = ~0 + infection_history, - truncate_lower = FALSE) + strict_lower_limit = FALSE) suppressWarnings({plot <- mod$plot_model_inputs()}) vdiffr::expect_doppelganger("inputdata_covariates", plot) }) diff --git a/vignettes/biokinetics.Rmd b/vignettes/biokinetics.Rmd index c97a926..9d735f5 100644 --- a/vignettes/biokinetics.Rmd +++ b/vignettes/biokinetics.Rmd @@ -28,7 +28,7 @@ The `fit` method then has the same function signature as the underlying [cmdstan ```{r, warning=FALSE, message=FALSE, results="hide"} dat <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) mod <- epikinetics::biokinetics$new(data = dat, - lower_detection_limit = 40, + lower_censoring_limit = 40, covariate_formula = ~0 + infection_history) delta <- mod$fit(parallel_chains = 4, iter_warmup = 50,