---
title: "Scenarios"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Scenarios}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

```{r, results='asis', echo=F, message=F, warning=F}
if (campsis::on_cran()) {
  cat(
    "This vignette was not built on CRAN. Please check out the online version [here](https://campsisverse.github.io/campsis.doc/articles/v11_scenarios.html)."
  )
  knitr::knit_exit()
}
```

```{r, results='hide', echo=F, message=F, warning=F}
library(campsis)
```

This vignette shows how scenarios can be implemented.

### Simulate several protocols at once

Using treatment arms:

```{r several_protocols_with_arms, fig.align='center', fig.height=4, fig.width=8, message=F}
model <- model_suite$nonmem$advan4_trans4

arm1 <- Arm(subjects = 3, label = "1000 mg SD") %>%
  add(Bolus(time = 0, 1000)) %>%
  add(Observations(times = seq(0, 24, by = 0.1)))

arm2 <- Arm(subjects = 5, label = "1500 mg SD") %>%
  add(Bolus(time = 0, 1500)) %>%
  add(Observations(times = seq(0, 24, by = 0.1)))

arm3 <- Arm(subjects = 10, label = "2000 mg SD") %>%
  add(Bolus(time = 0, 2000)) %>%
  add(Observations(times = seq(0, 24, by = 0.1)))

results <- simulate(
  model = model,
  dataset = Dataset() %>% add(c(arm1, arm2, arm3)),
  seed = 1
)
spaghetti_plot(results, "CONC") +
  ggplot2::facet_wrap(~ARM)
```

Using scenarios:

```{r several_protocols_with_scenarios, fig.align='center', fig.height=4, fig.width=8, message=F}
model <- model_suite$nonmem$advan4_trans4

dataset <- Dataset() %>%
  add(Observations(times = seq(0, 24, by = 0.1)))

scenarios <- Scenarios() %>%
  add(Scenario(
    "1000 mg SD",
    dataset = ~ .x %>% set_subjects(3) %>% add(Bolus(time = 0, 1000))
  )) %>%
  add(Scenario(
    "1500 mg SD",
    dataset = ~ .x %>% set_subjects(5) %>% add(Bolus(time = 0, 1500))
  )) %>%
  add(Scenario(
    "2000 mg SD",
    dataset = ~ .x %>% set_subjects(10) %>% add(Bolus(time = 0, 2000))
  ))

results <- simulate(
  model = model,
  dataset = dataset,
  scenarios = scenarios,
  seed = 1
)
spaghetti_plot(results, "CONC") + ggplot2::facet_wrap(~SCENARIO)
```

### Make a model parameter vary

Assume we want to test different values of `THETA_KA`:

```{r scenarios_make_theta_ka_vary, fig.align='center', fig.height=4, fig.width=8, message=F}
model <- model_suite$nonmem$advan4_trans4

ds <- Dataset(50) %>%
  add(Bolus(time = 0, amount = 1000)) %>%
  add(Observations(times = seq(0, 24, by = 0.1)))

scenarios <- Scenarios() %>%
  add(Scenario(
    "THETA_KA=1",
    model = ~ .x %>% replace(Theta(name = "KA", value = 1))
  )) %>%
  add(Scenario(
    "THETA_KA=3",
    model = ~ .x %>% replace(Theta(name = "KA", value = 3))
  )) %>%
  add(Scenario(
    "THETA_KA=6",
    model = ~ .x %>% replace(Theta(name = "KA", value = 6))
  ))

results <- model %>% simulate(dataset = ds, scenarios = scenarios, seed = 1)
shaded_plot(results, "CONC")
```

### Compare different distributions

Assume we want to compare different distributions of body weight `BW`:

```{r scenarios_compare_bw_distributions , fig.align='center', fig.height=6, fig.width=8, message=F}
model <- model_suite$nonmem$advan1_trans2 %>%
  replace(Equation("CL", "THETA_CL*exp(ETA_CL)*pow(BW/70, 0.75)")) %>%
  disable("IIV")

ds <- Dataset(50) %>%
  add(Bolus(time = 0, amount = 1000)) %>%
  add(Observations(times = seq(0, 24, by = 0.1))) %>%
  add(Covariate("BW", 70))

scenarios <- Scenarios() %>%
  add(Scenario("Constant BW")) %>%
  add(Scenario(
    "BW ∼ Uniform distribution",
    dataset = ~ .x %>%
      replace(Covariate("BW", UniformDistribution(min = 60, max = 80)))
  )) %>%
  add(Scenario(
    "BW ∼ Normal distribution",
    dataset = ~ .x %>%
      replace(Covariate("BW", NormalDistribution(mean = 70, sd = 10)))
  ))

results <- simulate(
  model = model,
  dataset = ds,
  scenarios = scenarios,
  seed = 1
)
shaded_plot(results, "CONC") +
  ggplot2::facet_wrap(~SCENARIO, ncol = 1)
```