-
Notifications
You must be signed in to change notification settings - Fork 0
/
time.qmd
178 lines (135 loc) · 7 KB
/
time.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
# On the nature of time {#sec-time}
```{r}
#| include: false
source("common.R")
```
There are two main sorts of treatment of time within odin models:
**Continuous time models**, as seen in the previous chapter, are described in terms of a set of [ordinary differential equations](https://en.wikipedia.org/wiki/Ordinary_differential_equation) (ODEs). At any point in time we can compute the rates of change of our quantity of interest, and using that and some initial condition, we can generate a time-series of the quantity or quantities involved.
**Discrete time models**, which will be the focus of the next section, are both more basic and more complicated depending on how you approach them. We assume that the system takes steps of some size `dt` and our equations will describe how each variable is updated -- how it changes from a value at time `t` to time `t + dt`. Sometimes this formulation is called a set of "recurrence equations".
This seems like a small difference but it is quite profound.
* Often the same system can be expressed as either a set of ODEs or as a set of recurrence equations but the two methods will behave differently (see below)
* It is not knowable which formulation will be "faster"; ODE solutions must take more steps where the solution is complex in order to smooth it out, but at the same time can take very large steps where the solution is simple
* ODEs must be **deterministic**; a unique smooth solution exists. Discrete time models are free to be stochastic, as it does not really matter what happens between `t` and `t + dt`
::: {.callout-note}
If you have used odin version 1, then discrete time models had a different (though fundamentally equivalent) formulation of time. See the [odin2 migration vignette](https://mrc-ide.github.io/odin2/articles/migrating.html) for more details.
:::
## A simple example where time matters {#sec-time-logistic}
```{r}
library(odin2)
library(dust2)
```
Consider just about the simplest interesting "system" of ODEs, the [logistic equation](https://en.wikipedia.org/wiki/Logistic_function#Logistic_differential_equation):
```{r}
logistic_ode <- odin({
deriv(x) <- r * x * (1 - x)
initial(x) <- x0
x0 <- parameter(0.01)
r <- parameter(0.1)
}, debug = TRUE)
logistic_ode
```
As in @sec-odin-sir, we can run this through time from its initial conditions:
```{r}
sys <- dust_system_create(logistic_ode, list())
dust_system_set_state_initial(sys)
t <- seq(0, 100, by = 1)
y <- dust_system_simulate(sys, t)
plot(t, y[1, ], type = "l", xlab = "Time", ylab = "Population")
```
Now, we write out a closely related system which uses discrete time; the "[logistic map](https://en.wikipedia.org/wiki/Logistic_map)"
```{r}
logistic_map <- odin({
update(x) <- r * x * (1 - x)
initial(x) <- x0
x0 <- parameter(0.01)
r <- parameter(1.1)
}, debug = TRUE)
```
and run it:
```{r}
sys2 <- dust_system_create(logistic_map, list())
dust_system_set_state_initial(sys2)
y2 <- dust_system_simulate(sys2, t)
plot(t, y[1, ], type = "l", xlab = "Time", ylab = "Population")
lines(t, y2[1, ], col = "red")
```
Let's write a function for this:
```{r}
run <- function(r, sys, t) {
dust_system_update_pars(sys, list(r = r))
dust_system_set_time(sys, 0)
dust_system_set_state_initial(sys)
dust_system_simulate(sys, t)[1, ]
}
```
We can run this function for some small values of `r` and see the population go extinct (black line) or grow to some stable value (red lines)
```{r}
plot(t, run(0.5, sys2, t), col = "black", ylim = 0:1, type = "l",
xlab = "Time", ylab = "Population")
lines(t, run(1.5, sys2, t), col = "red")
lines(t, run(1.8, sys2, t), col = "red")
lines(t, run(2, sys2, t), col = "red")
```
The equilibrium level of the discrete time system is not 0 or 1 (as in the continuous time version) but `(r - 1) / r`.
For larger values of `r` we start getting oscillations around the maximum, either decaying (blue lines) or stable (green line):
```{r}
plot(t, run(2.1, sys2, t), col = "blue", ylim = 0:1, type = "l",
xlab = "Time", ylab = "Population")
lines(t, run(2.5, sys2, t), col = "blue")
lines(t, run(2.9, sys2, t), col = "blue")
lines(t, run(3.3, sys2, t), col = "green")
```
Higher values bounce around chaotically:
```{r}
plot(t, run(3.8, sys2, t), ylim = 0:1, type = "l",
xlab = "Time", ylab = "Population")
```
None of these dynamics appeared in the continuous time version!
## Variables that reset {#sec-reset-variables}
Let's revisit our simple SIR example from @sec-odin-sir. Our model had just three compartments; `S`, `I` and `R`. These all represent prevalences, but frequently we are interested in incidences, for example the incidence of new infections or cases. This might be because we want to fit to data that are incidences, or it might just be that we want to run projections of a model and output incidences over time.
In some simple models it might be easy enough to extract incidences given outputs on all states of the model, but with increasingly complexity this calculation quickly becomes overly complicated or even impossible. Thus we support calculating incidences within your odin code.
Suppose that each time unit represents one day. We can declare `incidence` as a variable that resets every day by using the `zero_every` argument.
```{r}
sir <- odin({
deriv(S) <- -beta * S * I / N
deriv(I) <- beta * S * I / N - gamma * I
deriv(R) <- gamma * I
deriv(incidence) <- beta * S * I / N
initial(S) <- N - I0
initial(I) <- I0
initial(R) <- 0
initial(incidence, zero_every = 1) <- 0
N <- parameter(1000)
I0 <- parameter(10)
beta <- parameter(0.2)
gamma <- parameter(0.1)
})
```
If we were interested in weekly incidence we might choose to write
```r
initial(incidence, zero_every = 7) <- 0
```
or to work with time so that one unit represents a *week* rather than a day and scaling our rates accordingly.
::: {.callout-note}
Variables reset on multiples of the value of `zero_every`, and it is not currently possible to have this offset. If `zero_every = 7` the variable would reset only on multiples of 7 - so if you were dealing with weekly data, you would have to be careful to set time-zero in such a way that your weekly data aligns on multiples of 7!
:::
When running this model we'll see incidence produced:
```{r}
pars <- list(beta = 1, gamma = 0.6)
sys <- dust_system_create(sir, pars)
dust_system_set_state_initial(sys)
t <- seq(0, 20, by = 0.05)
y <- dust_system_simulate(sys, t)
y <- dust_unpack_state(sys, y)
plot(t, y$incidence, type = "s", xlab = "Time", ylab = "Incidence")
```
We have output at intervals of one 20th of a day, so we see 20 levels of incidence per day, with the last being the daily peak. If we output instead at intervals of one day, we only get those peaks corresponding to the daily incidence
```{r}
sys <- dust_system_create(sir, pars)
dust_system_set_state_initial(sys)
t <- seq(0, 20)
y <- dust_system_simulate(sys, t)
y <- dust_unpack_state(sys, y)
plot(t, y$incidence, type = "p", xlab = "Time", ylab = "Incidence")
```
We will see that resetting variables works in a similar fashion for discrete-time models in @sec-stochasticity.