# Get started

####

#### 2023-12-08

    library(dplyr)
    library(mrgsim.sa)
    library(patchwork)

mrgsim.sa is a package that allows you to do various types of ad-hoc or local sensitivity analysis for models written in mrgsolve. This vignette will help you get started.

By “ad-hoc” sensitivity analysis, I mean selecting some model parameters of interest, varying them at certain discrete, systematic way, simulating from those varied parameters and visualizing the outputs. For example, “what happens when you double the value of a parameter or cut it in half; lets plots outputs for parameters within those two extremes, filling in with 3 or 4 values in between.”

So first, we need a model. This can be any model you write using the mrgsolve package. Here, we’ll use the example `house` model provided by mrgsolve

    mod <- house(outvars = "CP,RESP")

This model has parameters

    param(mod)

    .
    .  Model parameters (N=14):
    .  name value . name  value
    .  CL   1     | SEX   0
    .  D1   2     | SEXCL 0.7
    .  F1   1     | SEXVC 0.85
    .  IC50 10    | VC    20
    .  KA   1.2   | WT    70
    .  KIN  100   | WTCL  0.75
    .  KOUT 2     | WTVC  1

and outputs

    outvars(mod)

    . $cmt
    . [1] "RESP"
    .
    . $capture
    . [1] "CP"

Let’s vary `CL` and `VC` and look at the `RESP` output.

As suggested in the example above, lets double and halve those parameters, and look a total of 5 values between those extremes. Use the `parseq_fct()` function after selecting the parameters you want to vary

    out <-
      mod %>%
      ev(amt = 100) %>%
      parseq_fct(CL, VC, .factor = 2, .n = 5) %>%
      sens_each()

The “base” value for each parameter is whatever is currently in the model; in this case it is

    param(mod)[c("CL", "VC")]

    . $CL
    . [1] 1
    .
    . $VC
    . [1] 20

Passing the `.factor` argument as 2 means multiply those base values by 2 for the upper extreme and 1/2 for the lower extreme. The `.n` argument says to fill in 5 parameter values between those two extremes.

The `sens_each()` function call above tells you to vary `CL` and `VC` one at a time.

The output is a tibble in long format, with class `sens_each`

    out

    . # A tibble: 9,680 × 7
    .     case  time p_name p_value dv_name dv_value ref_value
    .  * <int> <dbl> <chr>    <dbl> <chr>      <dbl>     <dbl>
    .  1     1  0    CL         0.5 RESP       50        50
    .  2     1  0    CL         0.5 RESP       50        50
    .  3     1  0    CL         0.5 CP          0         0
    .  4     1  0    CL         0.5 CP          0         0
    .  5     1  0    CL         0.5 RESP       50        50
    .  6     1  0    CL         0.5 RESP       50        50
    .  7     1  0    CL         0.5 CP          0         0
    .  8     1  0    CL         0.5 CP          0         0
    .  9     1  0.25 CL         0.5 RESP       48.7      48.7
    . 10     1  0.25 CL         0.5 CP          1.29      1.29
    . # ℹ 9,670 more rows

    class(out)

    . [1] "sens_each"  "tbl_df"     "tbl"        "data.frame"

Taking inventory of this output

    count(out, p_name, dv_name)

    . # A tibble: 4 × 3
    .   p_name dv_name     n
    .   <chr>  <chr>   <int>
    . 1 CL     CP       2420
    . 2 CL     RESP     2420
    . 3 VC     CP       2420
    . 4 VC     RESP     2420

We pass this output object to `sens_plot()` and name the variable we want to plot

    sens_plot(out, "RESP")

There are other ways to vary parameters in the model

- `parseq_fct()` - increase and decrease by a certain factor
- `parseq_cv()` - increase and decrease by a certain coefficient of variation
- `parseq_range()` - manually specify the range for varying parameters
- `parseq_manual()` - manually specify all values for the parameters

For example, to vary `CL` by 60% coefficient of variation, plotting 5 values between -2 and 2 sd and looking at `CP` output

    mod %>%
      ev(amt = 100) %>%
      parseq_cv(CL, .cv = 50, .nsd = 2) %>%
      sens_each() %>%
      sens_plot("CP")

Or we can look at how `VC` influences approach to steady state

    out <-
      mod %>%
      ev(amt = 100, ii = 24, addl = 10) %>%
      update(end = 240) %>%
      parseq_manual(VC = seq(10,100,20)) %>%
      sens_each()

    sens_plot(out, "CP")

We can also look at multiple outputs on the same plot

    out <-
      mod %>%
      ev(amt = 100) %>%
      parseq_fct(CL, KA, .n = 3) %>%
      sens_each(end = 36)

The `facet_wrap` layout puts parameters in rows and outputs in columns

    sens_plot(out, layout = "facet_wrap")

The `facet_grid` layout puts outputs in rows and parameters in columns

    sens_plot(out, layout = "facet_grid")

You can also plot this in “grid” format, where the actual parameter values are shown in the legend

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

Or look at multiple outputs

    out %>%
      select_sens(dv_name = "RESP,CP") %>%
      sens_plot(grid = TRUE) %>%
      patchwork::wrap_plots(ncol = 1)