diff --git a/DESCRIPTION b/DESCRIPTION index 81b40d8..4cc560f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -18,7 +18,9 @@ Imports: logger, mosaic, pkgload, + plotly, R6, + shiny, stats, stringr, tidybayes diff --git a/NAMESPACE b/NAMESPACE index 620d689..b500fc9 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -40,4 +40,6 @@ importFrom(ggplot2,scale_x_continuous) importFrom(ggplot2,scale_x_date) importFrom(ggplot2,scale_y_continuous) importFrom(ggplot2,sec_axis) +importFrom(ggplot2,theme) +importFrom(ggplot2,unit) useDynLib(epikinetics, .registration = TRUE) diff --git a/R/biokinetics.R b/R/biokinetics.R index 91487fb..885223f 100644 --- a/R/biokinetics.R +++ b/R/biokinetics.R @@ -347,37 +347,43 @@ biokinetics <- R6::R6Class( #' this plot is of the data as provided to the Stan model so is on a log scale, #' regardless of whether data was provided on a log or a natural scale. #' @param tmax Integer. Maximum time since last exposure to include. Default 150. + #' @param ncol Optional number of cols to display facets in. #' @return A ggplot2 object. - plot_model_inputs = function(tmax = 150) { + plot_model_inputs = function(tmax = 150, ncol = NULL) { plot_sero_data(private$data, tmax = tmax, + ncol = ncol, covariates = private$all_formula_vars, 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 Opens an RShiny app to help with model diagnostics. + inspect = function() { + inspect_model(self, private) + }, + #' @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() diff --git a/R/epikinetics-package.R b/R/epikinetics-package.R index 77fec4a..67e10e3 100644 --- a/R/epikinetics-package.R +++ b/R/epikinetics-package.R @@ -11,9 +11,9 @@ #' @importFrom data.table .NGRP #' @importFrom data.table .SD #' @importFrom data.table data.table -#' @importFrom ggplot2 aes annotate facet_wrap geom_point geom_ribbon geom_line geom_smooth geom_bar geom_density_2d -#' geom_vline geom_hline geom_path labs ggplot guides guide_legend scale_y_continuous -#' scale_x_continuous scale_x_date sec_axis +#' @importFrom ggplot2 aes annotate facet_wrap geom_point geom_ribbon geom_line geom_smooth geom_bar +#' geom_vline geom_hline geom_path labs ggplot guides guide_legend scale_y_continuous theme unit +#' geom_density_2d scale_x_continuous scale_x_date sec_axis #' @useDynLib epikinetics, .registration = TRUE ## usethis namespace: end diff --git a/R/inspect-model.R b/R/inspect-model.R new file mode 100644 index 0000000..75cea53 --- /dev/null +++ b/R/inspect-model.R @@ -0,0 +1,199 @@ +inspect_model <- function(mod, private) { + + prior_inputs <- function(name, description) { + mu <- paste("mu", name, sep = "_") + sigma <- paste("sigma", name, sep = "_") + shiny::div(shiny::fluidRow( + shiny::column(4, + description + ), + shiny::column(4, + shiny::fluidRow(class = "form-group", + shiny::tags$label(paste0("mean (", mu, ")"), class = "col-sm-6 col-form-label text-right"), + shiny::column(6, raw_numeric_input(mu, value = private$priors[[mu]])) + ) + ), + shiny::column(4, + shiny::fluidRow(class = "form-group", + shiny::tags$label(paste0("SD (", sigma, ")"), class = "col-sm-6 col-form-label text-right"), + shiny::column(6, raw_numeric_input(sigma, value = private$priors[[sigma]])), + ) + )) + ) + } + + all_covariates <- c("None", detect_covariates(private$data)) + + ui <- shiny::fluidPage(style = "margin: 0.5em", + shiny::fluidRow( + shiny::column(5, + shiny::h3("Prior predictive check"), + plotly::plotlyOutput("prior_predicted"), + shiny::tags$pre(style = "overflow: hidden; text-wrap: auto; word-break: keep-all; white-space: pre-line; margin-top: 20px;", + shiny::textOutput("prior_code", inline = TRUE) + ), + prior_inputs("t0", "Baseline titre value"), + prior_inputs("tp", "Time to peak titre"), + prior_inputs("ts", "Time to start of waning"), + prior_inputs("m1", "Boosting rate"), + prior_inputs("m2", "Plateau rate"), + prior_inputs("m3", "Waning rate") + ), + shiny::column(7, + shiny::h3("Model input data"), + shiny::uiOutput( + "data_plot" + ), + shiny::div(style = "margin-top: 20px;", + shiny::fluidRow(class = "form-group", + shiny::column(2, + shiny::numericInput("ncol", label = "Number of columns", value = 3) + ), + shiny::column(3, + shiny::selectInput("covariate", "Facet by", + choices = all_covariates, + selected = "None", + selectize = FALSE) + ), + shiny::column(7, + shiny::div(class = "form-group", + shiny:::shinyInputLabel("filter", "Filter by"), + shiny::fluidRow( + shiny::column(5, + raw_select_input("filter", + choices = all_covariates, + selected = "None") + ), + shiny::column(1, style = "padding-top: 5px;", "~="), + shiny::column(5, + raw_text_input("filter_value", placeholder = "regex") + ) + ) + ) + ) + ) + ) + ) + ), + shiny::fluidRow(style = "margin-top: 20px;", + shiny::column(12, + shiny::h3(shiny::textOutput("fitted")) + ) + ) + ) + + server <- function(input, output, session) { + # priors + prior <- shiny::reactive( + biokinetics_priors(mu_t0 = input$mu_t0, mu_tp = input$mu_tp, + mu_ts = input$mu_ts, mu_m1 = input$mu_m1, + mu_m2 = input$mu_m2, mu_m3 = input$mu_m3, + sigma_t0 = input$sigma_t0, sigma_tp = input$sigma_tp, + sigma_ts = input$sigma_ts, sigma_m1 = input$sigma_m1, + sigma_m2 = input$sigma_m2, sigma_m3 = input$sigma_m3) + ) + output$prior_code <- shiny::renderText({ + prior_code(input) + }) + output$prior_predicted <- plotly::renderPlotly({ + plotly::style(plotly::ggplotly(plot(prior(), + data = private$data, + upper_censoring_limit = private$stan_input_data$upper_censoring_limit, + lower_censoring_limit = private$stan_input_data$lower_censoring_limit)), textposition = "right") + }) + + # model inputs + cols <- shiny::reactive({ + if (is.na(input$ncol)) { + return(NULL) + } else { + return(input$ncol) + } + }) + + selected_covariate <- shiny::reactive({ + input$covariate + }) + + filter <- shiny::reactive({ + input$filter + }) + + filter_value <- shiny::reactive({ + input$filter_value + }) + + data <- shiny::reactive({ + if (filter_value() != "" && + !is.null(filter()) && + filter() != "None") { + return(private$data[grepl(filter_value(), get(filter()), ignore.case = TRUE)]) + } else { + return(private$data) + } + }) + + plot_inputs <- shiny::reactive({ + selected <- selected_covariate() + if (is.null(selected) || selected == "None") { + selected <- character(0) + } + plot_sero_data(data(), + ncol = cols(), + covariates = selected, + upper_censoring_limit = private$stan_input_data$upper_censoring_limit, + lower_censoring_limit = private$stan_input_data$lower_censoring_limit) + + theme(plot.margin = unit(c(1, 0, 0, 0), "cm")) + }) + + output$data <- plotly::renderPlotly({ + if (nrow(data()) > 0) { + gp <- plotly::style(plotly::ggplotly(plot_inputs()), textposition = "right") + if (selected_covariate() != "None") { + return(facet_strip_bigger(gp, 30)) + } else { + return(gp) + } + } + }) + + output$data_plot <- shiny::renderUI({ + if (nrow(data()) > 0) { + plotly::plotlyOutput("data") + } else { + shiny::h3("No rows selected. Please change your filter.") + } + }) + + # model outputs + output$fitted <- shiny::renderText({ + if (is.null(private$fitted)) { + "Model has not been fitted yet. Once fitted, inspect the model again to see posterior predictions." + } + }) + } + + logger::log_info( + "Starting Shiny app for model review; use Ctrl + C to quit" + ) + shiny::runApp( + shiny::shinyApp(ui, server), + quiet = TRUE, + launch.browser = shiny::paneViewer() + ) + invisible() +} + +# plotly can't handle multi-line facet titles, so manually make +# the facet titles a little bigger when there are covariates +facet_strip_bigger <- function(gp, size) { + + n_facets <- c(1:length(gp[["x"]][["layout"]][["shapes"]])) + + for (i in n_facets) { + gp[["x"]][["layout"]][["shapes"]][[i]][["y0"]] <- +as.numeric(size) + gp[["x"]][["layout"]][["shapes"]][[i]][["y1"]] <- 0 + } + + return(gp) +} diff --git a/R/plot.R b/R/plot.R index b15b595..83f59c9 100644 --- a/R/plot.R +++ b/R/plot.R @@ -56,10 +56,8 @@ plot.biokinetics_priors <- function(x, geom_point(data = dat, size = 0.5, aes(x = time_since_last_exp, y = value)) - - plot <- add_limits(plot, upper_censoring_limit, lower_censoring_limit) } - plot + add_limits(plot, upper_censoring_limit, lower_censoring_limit) } #' @title Plot serological data @@ -69,11 +67,13 @@ plot.biokinetics_priors <- function(x, #' @return A ggplot2 object. #' @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 ncol Integer. Optional number of columns to display facets in. #' @param covariates Optional vector of covariate names to facet by (these must correspond to columns in the data.table) #' @param upper_censoring_limit Optional upper detection limit. #' @param lower_censoring_limit Optional lower detection limit. plot_sero_data <- function(data, tmax = 150, + ncol = NULL, covariates = character(0), upper_censoring_limit = NULL, lower_censoring_limit = NULL) { @@ -87,7 +87,7 @@ plot_sero_data <- function(data, 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)))) + + facet_wrap(eval(parse(text = facet_formula(covariates))), ncol = ncol) + guides(colour = guide_legend(title = "Titre type")) add_limits(plot, upper_censoring_limit, lower_censoring_limit) @@ -305,7 +305,7 @@ add_limits <- function(plot, upper_censoring_limit, lower_censoring_limit) { linetype = 'dotted') + annotate("text", x = 1, y = lower_censoring_limit, - label = "Lower detection limit", + label = "Lower censoring limit", vjust = -0.5, hjust = 0, size = 3) @@ -316,7 +316,7 @@ add_limits <- function(plot, upper_censoring_limit, lower_censoring_limit) { linetype = 'dotted') + annotate("text", x = 1, y = upper_censoring_limit, - label = "Upper detection limit", + label = "Upper censoring limit", vjust = -0.5, hjust = 0, size = 3) diff --git a/R/shiny-utils.R b/R/shiny-utils.R new file mode 100644 index 0000000..dddba9f --- /dev/null +++ b/R/shiny-utils.R @@ -0,0 +1,44 @@ +raw_numeric_input <- function(inputId, value, min = NA, max = NA, step = NA) { + value <- shiny::restoreInput(id = inputId, default = value) + inputTag <- shiny::tags$input(id = inputId, type = "number", class = "shiny-input-number form-control", value = shiny:::formatNoSci(value)) + if (!is.na(min)) inputTag$attribs$min <- min + if (!is.na(max)) inputTag$attribs$max <- max + if (!is.na(step)) inputTag$attribs$step <- step + inputTag +} + +raw_text_input <- function(inputId, value = "", placeholder = NULL) { + value <- shiny::restoreInput(id = inputId, default = value) + shiny::tags$input(id = inputId, type = "text", class = "shiny-input-text form-control", value = value, placeholder = placeholder) +} + +raw_select_input <- function(inputId, choices, selected = NULL, multiple = FALSE) { + selected <- shiny::restoreInput(id = inputId, default = selected) + choices <- shiny:::choicesWithNames(choices) + if (is.null(selected)) { + if (!multiple) selected <- shiny:::firstChoice(choices) + } else selected <- as.character(selected) + shiny::tags$select(id = inputId, class = "shiny-input-select", class = "form-control", shiny:::selectOptions(choices, selected, inputId)) +} + + +prior_code <- function(input) { + deparse(substitute(biokinetics_priors(mu_t0 = a, mu_tp = b, + mu_ts = c, mu_m1 = d, + mu_m2 = e, mu_m3 = f, + sigma_t0 = g, sigma_tp = h, + sigma_ts = i, sigma_m1 = j, + sigma_m2 = k, sigma_m3 = l), + list(a = input$mu_t0, b = input$mu_tp, + c = input$mu_ts, d = input$mu_m1, + e = input$mu_m2, f = input$mu_m3, + g = input$sigma_t0, h = input$sigma_tp, + i = input$sigma_ts, j = input$sigma_m1, + k = input$sigma_m2, l = input$sigma_m3)), width.cutoff = 500L) +} + +detect_covariates <- function(data) { + setdiff(colnames(data), c("pid", "day", "last_exp_day", + "titre_type", "value", "censored", + "obs_id", "time_since_last_exp")) +} diff --git a/README.md b/README.md index f16d6ba..f3c28a1 100644 --- a/README.md +++ b/README.md @@ -81,6 +81,12 @@ To run all tests locally: devtools::test() ``` +To run tests in a single file: + +```{r} +devtools::test(filter="filename") +``` + Some tests are skipped on CI to avoid exorbitantly long build times, but this means it is important to run all tests locally at least once before merging a pull request. diff --git a/man/biokinetics.Rd b/man/biokinetics.Rd index b037694..e0db6fa 100644 --- a/man/biokinetics.Rd +++ b/man/biokinetics.Rd @@ -13,6 +13,7 @@ fit it to a dataset. \item \href{#method-biokinetics-new}{\code{biokinetics$new()}} \item \href{#method-biokinetics-plot_prior_predictive}{\code{biokinetics$plot_prior_predictive()}} \item \href{#method-biokinetics-plot_model_inputs}{\code{biokinetics$plot_model_inputs()}} +\item \href{#method-biokinetics-inspect}{\code{biokinetics$inspect()}} \item \href{#method-biokinetics-get_stan_data}{\code{biokinetics$get_stan_data()}} \item \href{#method-biokinetics-get_covariate_lookup_table}{\code{biokinetics$get_covariate_lookup_table()}} \item \href{#method-biokinetics-fit}{\code{biokinetics$fit()}} @@ -111,19 +112,31 @@ 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, regardless of whether data was provided on a log or a natural scale. \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{biokinetics$plot_model_inputs(tmax = 150)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{biokinetics$plot_model_inputs(tmax = 150, ncol = NULL)}\if{html}{\out{
}} } \subsection{Arguments}{ \if{html}{\out{
}} \describe{ \item{\code{tmax}}{Integer. Maximum time since last exposure to include. Default 150.} + +\item{\code{ncol}}{Optional number of cols to display facets in.} } \if{html}{\out{
}} } \subsection{Returns}{ A ggplot2 object. } +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-biokinetics-inspect}{}}} +\subsection{Method \code{inspect()}}{ +Opens an RShiny app to help with model diagnostics. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{biokinetics$inspect()}\if{html}{\out{
}} +} + } \if{html}{\out{
}} \if{html}{\out{}} diff --git a/man/plot_sero_data.Rd b/man/plot_sero_data.Rd index 6986bc0..ba4fcaa 100644 --- a/man/plot_sero_data.Rd +++ b/man/plot_sero_data.Rd @@ -7,6 +7,7 @@ plot_sero_data( data, tmax = 150, + ncol = NULL, covariates = character(0), upper_censoring_limit = NULL, lower_censoring_limit = NULL @@ -17,6 +18,8 @@ plot_sero_data( \item{tmax}{Integer. The number of time points in each simulated trajectory. Default 150.} +\item{ncol}{Integer. Optional number of columns to display facets in.} + \item{covariates}{Optional vector of covariate names to facet by (these must correspond to columns in the data.table)} \item{upper_censoring_limit}{Optional upper detection limit.} diff --git a/src/stan/antibody_kinetics_main.stan b/src/stan/antibody_kinetics_main.stan index 2d6cb38..06f6488 100644 --- a/src/stan/antibody_kinetics_main.stan +++ b/src/stan/antibody_kinetics_main.stan @@ -144,7 +144,7 @@ model { mu = boost_wane_ind( t, t0_ind, tp_ind, ts_ind, m1_ind, m2_ind, m3_ind, titre_type, id); - // Likelihood for uncensored observations + // Likelihood for uncensored observations (lower < value < upper) value[uncens_idx] ~ normal(mu[uncens_idx], sigma); // Likelihood for observations at lower limit diff --git a/tests/testthat/_snaps/plots/inputdata-covariates.svg b/tests/testthat/_snaps/plots/inputdata-covariates.svg index b13ff7f..0315399 100644 --- a/tests/testthat/_snaps/plots/inputdata-covariates.svg +++ b/tests/testthat/_snaps/plots/inputdata-covariates.svg @@ -503,9 +503,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -697,9 +697,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -894,9 +894,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -1389,9 +1389,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -1872,9 +1872,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -2072,9 +2072,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit diff --git a/tests/testthat/_snaps/plots/inputdata.svg b/tests/testthat/_snaps/plots/inputdata.svg index 5ed93b5..eebe3d4 100644 --- a/tests/testthat/_snaps/plots/inputdata.svg +++ b/tests/testthat/_snaps/plots/inputdata.svg @@ -681,9 +681,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -1339,9 +1339,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -2015,9 +2015,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit diff --git a/tests/testthat/_snaps/plots/populationtrajectories-data.svg b/tests/testthat/_snaps/plots/populationtrajectories-data.svg index 162ae5c..e1069aa 100644 --- a/tests/testthat/_snaps/plots/populationtrajectories-data.svg +++ b/tests/testthat/_snaps/plots/populationtrajectories-data.svg @@ -581,9 +581,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -799,9 +799,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -1023,9 +1023,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -1598,9 +1598,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -2152,9 +2152,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -2382,9 +2382,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit diff --git a/tests/testthat/_snaps/plots/populationtrajectories-logscale.svg b/tests/testthat/_snaps/plots/populationtrajectories-logscale.svg index 44e8c27..3b63f3d 100644 --- a/tests/testthat/_snaps/plots/populationtrajectories-logscale.svg +++ b/tests/testthat/_snaps/plots/populationtrajectories-logscale.svg @@ -581,9 +581,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -799,9 +799,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -1023,9 +1023,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -1598,9 +1598,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -2152,9 +2152,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit @@ -2382,9 +2382,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit diff --git a/tests/testthat/_snaps/plots/priorpredictive.svg b/tests/testthat/_snaps/plots/priorpredictive.svg index bef47bb..7246078 100644 --- a/tests/testthat/_snaps/plots/priorpredictive.svg +++ b/tests/testthat/_snaps/plots/priorpredictive.svg @@ -2282,9 +2282,9 @@ -Lower detection limit +Lower censoring limit -Upper detection limit +Upper censoring limit diff --git a/tests/testthat/test-convert-log-scale.R b/tests/testthat/test-convert-log-scale.R index 76ad408..3b88db9 100644 --- a/tests/testthat/test-convert-log-scale.R +++ b/tests/testthat/test-convert-log-scale.R @@ -2,7 +2,6 @@ test_that("Can convert to and from log scale in R", { inputs <- data.table::fread(test_path("testdata", "testdata.csv")) log_inputs <- convert_log2_scale(inputs, "me", smallest_value = 2) unlog_inputs <- convert_log2_scale_inverse(log_inputs, "me", smallest_value = 2) - expect_equal(inputs, unlog_inputs) }) diff --git a/tests/testthat/test-priors.R b/tests/testthat/test-priors.R index 1ca890e..9997fdb 100644 --- a/tests/testthat/test-priors.R +++ b/tests/testthat/test-priors.R @@ -1,5 +1,5 @@ -test_that("Can construct cab prior parameters", { - priors <- biokinetics_priors(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12) +test_that("Can construct biokinetics priors", { + priors <- biokinetics_priors(1, 2, 3, 4, 5, 6 , 7, 8, 9, 10, 11, 12) expect_s3_class(priors, "biokinetics_priors") expect_true(is.list(priors)) expect_equal(unclass(priors), list("mu_t0" = 1, "mu_tp" = 2, "mu_ts" = 3, diff --git a/tests/testthat/test-shiny-utils.R b/tests/testthat/test-shiny-utils.R new file mode 100644 index 0000000..23b1934 --- /dev/null +++ b/tests/testthat/test-shiny-utils.R @@ -0,0 +1,6 @@ +test_that("Can create code snippet with priors", { + expect_equal(prior_code(biokinetics_priors()), + paste("biokinetics_priors(mu_t0 = 4, mu_tp = 10, mu_ts = 60, mu_m1 = 0.25,", + "mu_m2 = -0.02, mu_m3 = 0, sigma_t0 = 2, sigma_tp = 2, sigma_ts = 3,", + "sigma_m1 = 0.01, sigma_m2 = 0.01, sigma_m3 = 0.01)")) +}) diff --git a/vignettes/diagnostics.Rmd b/vignettes/diagnostics.Rmd new file mode 100644 index 0000000..3268e4f --- /dev/null +++ b/vignettes/diagnostics.Rmd @@ -0,0 +1,64 @@ +--- +title: "Model diagnostics and plotting" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Model diagnostics and plotting} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +## Plotting model inputs + +### Priors +Calling `plot` on an object of type [biokinetics_priors](../reference/biokinetics_priors.html) will generate +a plot of the kinetics predicted by the given priors. + +```{r} +priors <- epikinetics::biokinetics_priors() +plot(priors) +``` + +You can optionally pass a dataset to compare against the predicted kinetics. Required columns are +`time_since_last_exp` and `value` and values should be on the scale required by the model. + +```{r} +data <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) +data[, time_since_last_exp := as.integer(day - last_exp_day, units = "days")] +data <- epikinetics::convert_log2_scale(data, min(data$value)) +priors <- epikinetics::biokinetics_priors() +plot(priors, data = data) +``` + +If you have an instance of the [biokinetics](../reference/biokinetics.html) class, the method `plot_prior_predictive` +generates this plot for the priors and data given to the model. + +```{r} +data <- data.table::fread(system.file("delta_full.rds", package = "epikinetics")) +priors <- epikinetics::biokinetics_priors() +mod <- epikinetics::biokinetics$new(priors = priors, data = data) +mod$plot_prior_predictive() +``` + +### Data + +If you have an instance of the [biokinetics](../reference/biokinetics.html) class, the method `plot_model_inputs` +plots the input data used to fit the model, disaggregated by the covariates in the covariate formula. + +```{r} +mod <- epikinetics::biokinetics$new(priors = priors, data = data, covariate_formula = ~0 + infection_history) +mod$plot_model_inputs() +``` + +## Interactive data exploration + +To play around with different priors and visualise input data filtered and disaggregated in different ways, +the function [biokinetics$inspect](reference/biokinetics.html#method-biokinetics-inspect) runs a local RShiny app with interactive plots. + +![RShiny demonstration](./shiny.webm) diff --git a/vignettes/shiny.webm b/vignettes/shiny.webm new file mode 100644 index 0000000..75e2122 Binary files /dev/null and b/vignettes/shiny.webm differ