part3.qmd

---
title: "Part 3: Multiple Runs"
---

## Introduction

The purpose of the `run_multiple` function is to run a specified number of simulations using the same model object. That is, this function makes it possible to compare model results across several separate and repeated simulations.

## The Principle Behind Multiple Runs

```{r, echo = FALSE}
# Set the seed for reproducibility
set.seed(123)

# Define the range of the uniform distribution
lower <- 0
upper <- 6

# Initialize vectors to store the sample sizes and their corresponding averages
sample_sizes <- c()
averages <- c()

# Set the maximum sample size
max_n <- 1000

# Generate samples and calculate averages for increasing sample sizes
for (n in 1:max_n) {
  # Generate a random sample from the uniform distribution
  sample <- runif(n, min = lower, max = upper)
  # Calculate the average of the sample
  average <- mean(sample)
  # Store the sample size and the corresponding average
  sample_sizes <- c(sample_sizes, n)
  averages <- c(averages, average)
}
plot(sample_sizes, averages, type = "l", xlab = "Sample Size (n)", 
     ylab = "Average", main = "Average vs. Sample Size")
```

In statistics, the Law of Large Numbers ensures that as the sample size increases, the sample mean (average) of a random variable will converge to the population mean. The same principle applies when using multiple runs to simulate epidemiological models. As the number of `epiworldR` simulations increases, the sample means of the reproductive number or model parameters, for example, will converge to their corresponding population means.

## Example: Simulating a SEIRCONN Model 500 Times

#### Setup and Running Model

To use the `run_multiple` function in `epiworldR`, create your `epiworldR_model` of choice; in this case, the example uses a `SEIRCONN` model for Measles, 10000 people, an initial prevalence of 0.0001 (0.01%), a contact rate of 2, probability of transmission 0.5, a total of 7 incubation days, and probability of recovery $\frac{1}{7}$.

```{r sirconn-setup}
library(epiworldR)
model_seirconn <- ModelSEIRCONN(
  name              = "Measles",
  n                 = 5000, 
  prevalence        = 0.001, 
  contact_rate      = 2,
  transmission_rate = 0.5,
  incubation_days   = 7,
  recovery_rate     = 1/7
  )
```

#### Generating a Saver

Next, generate a saver for the purpose of extracting the `total_hist` and `reproductive` information from the model object. Keep in mind that you can generate a saver for any metric compatible with the `make_saver` function (see details section of the `make_saver` help manual).

```{r}
# Generating a saver
# ?make_saver
saver <- make_saver("total_hist", "reproductive")
```

#### Running the Simulation

Now, use the `run_multiple` function with the model object, number of desired days to run the simulation, number of simulations to run, the generated saver, and number of threads for parallel computing.

```{r saver-generation}
# Running and printing
run_multiple(model_seirconn, ndays = 50, nsims = 500, saver = saver, nthread = 4)
```

#### Extracting Results

Using the `run_multiple_get_results` function, extract the results from the model object that was simulated 500 times for comparison across simulations.

```{r extracting-results}
head(run_multiple_get_results(model_seirconn)$total_hist)
head(run_multiple_get_results(model_seirconn)$reproductive)
```

#### Plotting

To plot the model parameters and reproductive numbers over time using boxplots, extract the results from the model object using `run_multiple_get_results`. For this example, the dates are filtered to observe the model parameters over the first 20 days. Notice each boxplot in the below table represents the observed values from each of the 500 simulations for each date.

```{r plotting-seirconn-epicurves}
seirconn_500 <- run_multiple_get_results(model_seirconn)$total_hist
seirconn_500 <- seirconn_500[seirconn_500$date <= 20,]
plot(seirconn_500)
```

To view the a plot of the reproductive number over all 50 days for each of the 500 simulations, store the reproductive results to a new object using `run_multiple_get_results`, then plot using the `plot` function. Notice each source exposure date displays a boxplot representing the distribution of reproductive numbers across all 500 simulations for each date. As expected, the reproductive number on average, decreases over time.

```{r reproductive-number-plot}
seirconn_500_r <- run_multiple_get_results(model_seirconn)$reproductive
plot(seirconn_500_r)
```

## Exercise 1

Consider for this exercise that there is a Hepatitis A outbreak. Your goal is to observe the average reproductive number over 100 simulations. Using a `run_multiple` simulation, what is the average reproductive number over the course of the first 20 days? Use a `SEIRCONN` model with:\
- n = 10000\
- prevalence = 0.01\
- contact_rate = 2\
- transmission_rate = 0.5\
- incubation_days = 2\
- recovery_rate = $\frac{1}{7}$

::: {.callout-tip collapse="true"}
General Steps:\
1. Create `epiworldR_model`\
2. Generate saver\
3. Run `epiworldR_model`\
4. Plot average reproductive number
:::

```{r}
# Your solution here

```

## Exercise 2

Using the same `SEIRCONN` model from exercise 1, simulate a vaccine intervention for the previous exercise's Hepatitis A virus outbreak where 50% of individuals in the population will receive the vaccine on day 10. How then, does the average reproductive number behave over 20 days and 100 simulations? Assume the following parameters:\
- susceptibility_reduction = .9\
- transmission_reduction = .5\
- recovery_enhancer = .5\
- death_reduction = .9

::: {.callout-tip collapse="true"}
General Steps:\
1. Create tool\
2. Use `globalaction_tool` & `add_global_action`\
3. Generate saver & `run_multiple`\
4. Plot average reproductive number
:::

```{r}
# Your solution here

```