In kylebaron/mrgsim.sa: Sensitivity Analysis with 'mrgsolve'

A simple, clean workflow for sensitivity analysis with mrgsolve.

knitr::opts_chunk$set(
  comment = ".",
  fig.path = "man/figures/README-", 
  fig.width = 8.5
)
options(tibble.print_min = 5, tibble.print_max = 5)

library(mrgsim.sa)

mod <- mread("pk1", modlib(), end = 48, delta = 0.1)

param(mod)

PK model sensitivity analysis by factor

The nominal (in model) parameter value is divided and multiplied by a factor, generating minimum and maximum bounds for simulating a sequence of parameter values

out <- 
  mod %>% 
  ev(amt = 100) %>% 
  select_par(CL, V) %>% 
  parseq_fct(.n=8) %>% 
  sens_each() 

sens_plot(out, "CP")

The simulated data is returned in a long format

out

And you can plot with more informative color scale and legend

sens_plot(out, "CP", grid = TRUE)

HIV viral dynamic model

We look at latent infected cell pool development over ten years at different "burst" size, or the number of HIV particles released when one cell lyses.

mod <- mread("hiv", "inst/example")

mod %>% 
  update(end = 365*10) %>%
  parseq_range(N = c(900,1500), .n = 10) %>%
  sens_each(tscale = 1/365) %>% 
  sens_plot("L", grid = TRUE)

Sensitivity analysis on custom sequences

The model is rifampicin PBPK.

mod <- mread("inst/example/rifampicin.cpp") %>% update(delta = 0.1)

mod %>% 
  ev(amt = 600) %>% 
  parseq_manual(
    SFKp = seq_fct(.$SFKp, n = 20), 
    Kp_muscle = seq_even(0.001, 0.1, n = 6)
  ) %>% 
  sens_each() %>% 
  sens_plot("Ccentral")

Simulate a grid

To this point, we have always used sens_each so that each value for each parameter is simulated one at a time. Now, simulate the grid or all combinations.

We use parseq_cv here, which generates lower and upper bounds for the range using 50% coefficient of variation.

out <- 
  mod %>% 
  update(outvars = "Ccentral") %>%
  ev(amt = 600) %>% 
  parseq_cv(fBCLint_all_kg, .n = 7) %>% 
  parseq_cv(SFKp, Kp_muscle, .n = 3) %>% 
  sens_grid(recsort = 3) 

out

out %>% sens_plot("Ccentral")

Local sensitivity analysis

mod <- modlib("pk2", delta = 0.1, end = 72)

doses <- ev(amt = 100)

out <- lsa(mod, var = "CP", par = "CL,V2,Q", events = doses)

out

lsa_plot(out, pal = NULL)