Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .Rbuildignore
Original file line number Diff line number Diff line change
Expand Up @@ -13,3 +13,4 @@ Makefile
inst/covr
^\.github$
^_pkgdown\.yml$
CLAUDE\.md
72 changes: 72 additions & 0 deletions CLAUDE.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,72 @@
# CLAUDE.md

This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.

## Commands

```bash
# Run tests
make test
# or: Rscript -e "testthat::test_dir('tests/testthat')"

# Run a single test file
Rscript -e "testthat::test_file('tests/testthat/test-sens_each.R')"

# Generate documentation (roxygen2)
make doc
# or: Rscript -e 'devtools::document(".")'

# Build package (no vignettes)
make build

# Install from built tarball
make install

# Full check (no vignettes)
make check

# Code coverage
make covr

# Render README
make readme
```

## Architecture

This is an R package that adds sensitivity analysis workflows on top of [mrgsolve](https://mrgsolve.org/) ODE models. The core idea: attach parameter sequences to a model object, then simulate while varying those parameters.

### Source file load order (from DESCRIPTION Collate)

`utils.R` → `parseq.R` → `sens.R` → `AAA.R` → `lsa.R` → `sens_each.R` → `sens_grid.R` → `sens_plot.R` → `sens_run.R` → `seq.R`

### Key concepts

**Parameter sequences** (`parseq.R`, `seq.R`): Functions like `parseq_fct()`, `parseq_cv()`, `parseq_range()`, `parseq_manual()` attach a `sens_values` attribute to the mrgsolve model object. `seq_geo()`, `seq_fct()`, `seq_even()`, `seq_cv()` generate the underlying numeric sequences.

**Sensitivity simulation** (`sens_each.R`, `sens_grid.R`): `sens_each()` varies one parameter at a time (OAT); `sens_grid()` does full factorial combinations via `expand.grid()`. Both convert parameter lists to `idata` (individual-level data) and call mrgsolve's simulation engine. Output is a nested tibble with columns `p_name`, `.value`, and `data`.

**Visualization** (`sens_plot.R`): `sens_plot()` dispatches on `sens_each` or `sens_grid` class. Supports layouts: default (patchwork grid), `facet_grid`, `facet_wrap`, `list`. Reference case (nominal parameter values) plotted as black dashed line.

**Local sensitivity analysis** (`lsa.R`): `lsa()` perturbs each parameter by `eps` (default 1e-7) and computes finite-difference sensitivities. Output has its own `plot.lsa()` S3 method.

**Convenience wrapper** (`sens_run.R`): `sens_run()` accepts `method=` ("factor", "cv", "range", "manual") and `vary=` ("each", "grid") to combine parseq + sens in one call.

### Typical workflow

```r
mod %>%
select_par(CL, VC) %>% # tidyselect on model parameters
parseq_fct(.n = 8) %>% # attach sequences (factor method)
sens_each() %>% # OAT simulation → nested tibble
sens_plot("CP") # visualize dependent variable
```

### S3 classes

- `sens_each` / `sens_grid` / `sens_data`: output of simulation functions; have `as.data.frame()`, `as_tibble()`, `print()`, and `sens_plot()` methods
- `lsa`: output of `lsa()`; has `plot()` method

### Example models

`inst/example/` contains three mrgsolve `.cpp` model files (`hiv.cpp`, `gcsf.cpp`, `rifampicin.cpp`) used in tests and documentation.
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ Encoding: UTF-8
Language: en-US
Depends: mrgsolve
Roxygen: list(markdown = TRUE)
RoxygenNote: 7.2.3
RoxygenNote: 7.3.3
VignetteBuilder: knitr
Collate:
'utils.R'
Expand Down
1 change: 1 addition & 0 deletions Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -78,5 +78,6 @@ clean:
rm -f README.html
rm -rf DOCS
rm -rf mrgsim.sa.Rcheck
rm -rf mrgsim.sa*.900*.tar.gz


2 changes: 1 addition & 1 deletion R/AAA.R
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
#' @importFrom ggplot2 scale_color_viridis_c scale_color_brewer
#' @importFrom patchwork wrap_plots
#' @importFrom rlang quos sym set_names enexpr := as_string .env abort warn
#' @importFrom rlang abort is_integerish is_named
#' @importFrom rlang abort is_integerish is_named
#' @importFrom tidyr unnest nest pivot_longer
#' @importFrom tibble tibble as_tibble
#' @importMethodsFrom mrgsolve as.list param update as.data.frame
Expand Down
59 changes: 47 additions & 12 deletions R/sens_plot.R
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ sens_plot.sens_each <- function(data, dv_name = NULL, p_name = NULL,
logy = FALSE,
ncol = NULL, lwd = 0.8,
digits = 3, plot_ref = TRUE,
xlab = "time", ylab = dv_name[1],
xlab = "time", ylab = NULL,
layout = c("default", "facet_grid",
"facet_wrap", "list"),
grid = FALSE, ...) {
Expand All @@ -104,6 +104,20 @@ sens_plot.sens_each <- function(data, dv_name = NULL, p_name = NULL,
dv_name <- cvec_cs(dv_name)
}

if(is.null(ylab)) {
ylab <- dv_name
}

assert_that(is.character(xlab))
xlab <- xlab[1]

if(length(dv_name) != length(ylab)) {
ndv <- length(dv_name)
ny <- length(ylab)
msg <- glue("`dv_name` ({ndv}) and `ylab` ({ny}) have different lengths.")
abort(glue(msg))
}

if((grid && length(dv_name) > 1)) {
list <- TRUE
}
Expand All @@ -113,7 +127,7 @@ sens_plot.sens_each <- function(data, dv_name = NULL, p_name = NULL,
if(list) {
args <- c(as.list(environment()), list(...))
args$layout <- "default"
return(sens_plot_list(dv_name, args))
return(sens_plot_list(dv_name, ylab, args))
}

if(!is.null(p_name)) {
Expand Down Expand Up @@ -237,23 +251,28 @@ sens_plot.sens_each <- function(data, dv_name = NULL, p_name = NULL,
return(plots)
}

sens_plot_list <- function(dv_name, args) {
args$ylab <- NULL
sens_plot_list <- function(dv_name, ylab, args) {
out <- vector(mode = "list", length = length(dv_name))
i <- 1
for(this_dv_name in dv_name) {
args$dv_name <- this_dv_name
for(i in seq_along(dv_name)) {
args$dv_name <- dv_name[i]
args$ylab <- ylab[i]
out[[i]] <- do.call(sens_plot.sens_each, args)
i <- i+1
}
return(out)
}

#' @rdname sens_plot
#' @export
sens_plot.sens_grid <- function(data, dv_name = NULL, digits = 2, ncol = NULL,
lwd = 0.8, logy = FALSE,
plot_ref = TRUE, ...) { #nocov start
sens_plot.sens_grid <- function(data,
dv_name = NULL,
logy = FALSE,
ncol = NULL,
lwd = 0.8,
digits = 2,
plot_ref = TRUE,
xlab = "time",
ylab = dv_name,
...) { #nocov start

if(is.null(dv_name)) {
dv_name <- unique(data[["dv_name"]])
Expand All @@ -262,10 +281,25 @@ sens_plot.sens_grid <- function(data, dv_name = NULL, digits = 2, ncol = NULL,
dv_name <- cvec_cs(dv_name)
}

if(is.null(ylab)) {
ylab <- dv_name
}

assert_that(is.character(xlab))
xlab <- xlab[1]

if(length(dv_name) != length(ylab)) {
ndv <- length(dv_name)
ny <- length(ylab)
msg <- glue("`dv_name` ({ndv}) and `ylab` ({ny}) have different lengths.")
abort(glue(msg))
}

if(length(dv_name) > 1) {
args <- c(as.list(environment()), list(...))
out <- lapply(dv_name, function(this_dv_name) {
out <- Map(dv_name, ylab, f = function(this_dv_name, this_ylab) {
args$dv_name <- this_dv_name
args$ylab <- this_ylab
do.call(sens_plot.sens_grid, args)
})
return(out)
Expand Down Expand Up @@ -305,6 +339,7 @@ sens_plot.sens_grid <- function(data, dv_name = NULL, digits = 2, ncol = NULL,
p <- ggplot(data = data, aes(!!x, !!y, group=!!group, col=factor(!!group)))
p <- p + geom_line(lwd=lwd) + scale_color_discrete(name = pars[1])
p <- p + theme_bw() + theme(legend.position = "top")
p <- p + xlab(xlab) + ylab(ylab)
if(npar==2) p <- p + facet_wrap(formula, ncol = ncol)
if(npar==3) p <- p + facet_grid(formula)
if(isTRUE(logy)) p <- p + scale_y_log10()
Expand Down
8 changes: 5 additions & 3 deletions man/sens_plot.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

69 changes: 69 additions & 0 deletions tests/testthat/test-sens-plot.R
Original file line number Diff line number Diff line change
Expand Up @@ -111,3 +111,72 @@ test_that("sens_grid - multiple plots everything", {
expect_length(p, length(outv))
expect_is(p[[1]], "gg")
})

# xlab / ylab - sens_each -----------------------------------------------

test_that("sens_each - xlab is applied", {
p <- sens_plot(s1, "CP", xlab = "my-xlab")
expect_equal(p$labels$x, "my-xlab")
})

test_that("sens_each - ylab is applied for single dv", {
p <- sens_plot(s1, "CP", ylab = "my-ylab")
expect_equal(p$labels$y, "my-ylab")
})

test_that("sens_each - ylab defaults to dv_name for single dv", {
p <- sens_plot(s1, "CP")
expect_equal(p$labels$y, "CP")
})

test_that("sens_each - vectorized ylab applied to each plot in list", {
p <- sens_plot(s1, dv_name = "GUT,CP,RESP", ylab = c("gut", "conc", "resp"))
expect_is(p, "list")
expect_length(p, 3)
expect_equal(p[[1]]$labels$y, "gut")
expect_equal(p[[2]]$labels$y, "conc")
expect_equal(p[[3]]$labels$y, "resp")
})

test_that("sens_each - ylab length mismatch with dv_name errors", {
expect_error(
sens_plot(s1, dv_name = "GUT,CP,RESP", ylab = c("a", "b")),
regexp = "dv_name.*ylab.*different lengths"
)
})

test_that("sens_each - xlab must be character", {
expect_error(sens_plot(s1, "CP", xlab = 123), "xlab is not a character")
})

# xlab / ylab - sens_grid -----------------------------------------------

test_that("sens_grid - xlab is applied", {
p <- sens_plot(s2, "CP", xlab = "my-xlab")
expect_equal(p$labels$x, "my-xlab")
})

test_that("sens_grid - ylab is applied for single dv", {
p <- sens_plot(s2, "CP", ylab = "my-ylab")
expect_equal(p$labels$y, "my-ylab")
})

test_that("sens_grid - ylab defaults to dv_name for single dv", {
p <- sens_plot(s2, "CP")
expect_equal(p$labels$y, "CP")
})

test_that("sens_grid - vectorized ylab applied to each plot in list", {
p <- sens_plot(s2, dv_name = "CP,RESP", ylab = c("conc", "resp"))
expect_is(p, "list")
expect_length(p, 2)
expect_equal(p[[1]]$labels$y, "conc")
expect_equal(p[[2]]$labels$y, "resp")
})

test_that("sens_grid - ylab length mismatch with dv_name errors", {
expect_error(
sens_plot(s2, dv_name = "CP,RESP", ylab = "only-one"),
regexp = "dv_name.*ylab.*different lengths"
)
})