# 4 Simulation

All previous chapters were concerned with defining a compartmental model. In this chapter we switch to getting results from a defined model.

Once a model object is defined, it can be used to generate simulations using the `simulation_history`

function.

The output contains a column for the simulation date, one column for each state variable (`S`

, `I`

and `R`

in this case), and one column for every time-varying rate (`S_to_I`

). The names of the time-varying rates are always of the form `{from_state}_to_{to_state}`

. The reason why `S_to_I`

is time-varying in this model is that it depends on a state variable, `I`

, which is itself varying at every time-step. The `rate_summary`

function can be used to remind us of this fact.

```
## from to formula
## S_to_I S I (I) * (beta) * (1/N)
## I_to_R I R (gamma)
```

We see here that `S_to_I`

does indeed depend on `I`

in its formula, whereas `I_to_R`

depends only on a parameter, `gamma`

.

Note that the above command uses the tidyverse-style pipe, `%>%`

, operator and another tidyverse function, `select`

. This illustrates a general philosophy of McMasterPandemic, which is that we try to make the outputs plug into other existing and popular tools rather than reinvent existing functionality for a narrower purpose. For example, the `rate_summary`

function returns a data frame that can be manipulated by other data frame manipulation tools.

We can plug into other existing and popular tools to make a plot of the simulated epidemic trajectory.

```
(sir
%>% simulation_history
%>% select(-S_to_I)
%>% pivot_longer(-Date, names_to = "State", values_to = "Population")
%>% mutate(State = factor(State, levels = topological_sort(sir)))
%>% ggplot
+ geom_line(aes(Date, Population, colour = State))
)
```

There are a few places you can go from here:

- Learn how to fit a model to observed data through Calibration
- Learn how to modify the values of parameters in simulation time using Time Varying Parameters
- Keep reading to learn about simulating with Observation Error

## 4.1 Observation Error

```
sir_with_obs_err = (sir
%>% update_params(c(
nb_disp_S = 1e4,
nb_disp_I = 1e4,
nb_disp_R = 1e4
))
%>% update_error_dist(
S ~ negative_binomial("nb_disp_S"),
I ~ negative_binomial("nb_disp_I"),
R ~ negative_binomial("nb_disp_R")
)
)
```

```
set.seed(1L)
(sir_with_obs_err
%>% simulation_history(obs_error = TRUE)
%>% select(Date, S, I, R)
%>% pivot_longer(-Date, names_to = "var", values_to = "value")
%>% rename(date = Date)
%>% mutate(var = factor(var, levels = topological_sort(sir)))
%>% ggplot
+ facet_wrap(~ var)
+ geom_point(aes(date, value))
)
```