Skip to content

Commit b79f990

Browse files
stevencarlislewalkerstevencarlislewalker
committed
[skip ci] si with behavioural protection
1 parent 7cedf77 commit b79f990

File tree

7 files changed

+361
-0
lines changed

7 files changed

+361
-0
lines changed

‎inst/starter_models/si_behaviour/README.Rmd‎

Lines changed: 163 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,163 @@
1+
---
2+
title: "SI with Behaviour"
3+
index_entry: "SI model with behavioural protection in response to the number of infectious individuals"
4+
bibliography: ../../references.bib
5+
link-citations: TRUE
6+
author: Steve Walker
7+
output:
8+
github_document:
9+
toc: true
10+
editor_options:
11+
chunk_output_type: console
12+
---
13+
14+
```{r, include = FALSE}
15+
knitr::opts_chunk$set(
16+
collapse = TRUE,
17+
comment = "#>",
18+
fig.path = "./figures/"
19+
)
20+
theme_bw = function() ggplot2::theme_bw(base_size = 14)
21+
```
22+
23+
24+
This is an SI model in which individuals can protect themselves.
25+
26+
```{r flow_diagram, echo = FALSE, fig.height=3, fig.width=5, message=FALSE, warning=FALSE}
27+
library(macpan2)
28+
library(ggplot2)
29+
system.file("utils", "box-drawing.R", package = "macpan2") |> source()
30+
si = mp_tmb_library("starter_models", "si_behaviour", package = "macpan2")
31+
(si
32+
|> mp_layout_paths()
33+
|> plot_flow_diagram(show_flow_rates = TRUE)
34+
)
35+
```
36+
37+
The per-capita protection-rate, `alpha_t`, is time-varying, and increases as the number of infectious individuals, `I`, increases.
38+
39+
The code in this article uses the following packages.
40+
41+
```{r packages, message=FALSE, warning=FALSE}
42+
library(macpan2)
43+
library(ggplot2)
44+
library(dplyr)
45+
library(tidyr)
46+
```
47+
48+
49+
This model was inspired by a question from [Irena Papst](https://github.com/papsti). The [awareness model](https://github.com/canmod/macpan2/tree/main/inst/starter_models/awareness) provides another example of a behavioural response.
50+
51+
# States
52+
53+
54+
| Variable | Description |
55+
|------------|---------------------------------------------------------|
56+
| \( S(t) \) | Number of susceptible individuals |
57+
| \( I(t) \) | Number of infectious individuals |
58+
| \( P(t) \) | Number of individuals who have protected themselves |
59+
60+
61+
62+
# Parameters
63+
64+
| Parameter | Description |
65+
|------------------------|--------------------------------------------------------------------------|
66+
| \( \beta \) | Transmission rate |
67+
| \( \alpha_{\text{max}} \) | Maximum per-capita protection rate |
68+
| \( I^* \) | Threshold number of infectious individuals that triggers protection |
69+
| \( k \) | Steepness of the sigmoidal switching function |
70+
| \( N \) | Total population size (constant) |
71+
72+
73+
# Dynamics
74+
75+
$$
76+
\begin{align*}
77+
\sigma(I) &= \frac{1}{1 + \exp\left( -k \cdot (I - I^*) \right)} \\
78+
\alpha(t) &= \alpha_{\text{max}} \cdot \sigma(I(t)) \\
79+
\\
80+
\frac{dS}{dt} &= - \beta \cdot \frac{I}{N} \cdot S - \alpha(t) \cdot S \\
81+
\frac{dI}{dt} &= \beta \cdot \frac{I}{N} \cdot S \\
82+
\frac{dP}{dt} &= \alpha(t) \cdot S
83+
\end{align*}
84+
$$
85+
86+
87+
# Model Specification
88+
89+
This model has been specified in the `si_behaviour` directory [here](https://github.com/canmod/macpan2/blob/main/inst/starter_models/si_behaviour/tmb.R) and is accessible from the `macpan2` model library (see [Example Models](https://canmod.github.io/macpan2/articles/example_models.html) for details).
90+
91+
Printing the steps of the simulation loop illustrates how to implement a sigmoidal response of the protection rate to the prevalence of infectious individuals.
92+
93+
```{r flow_steps}
94+
mp_print_during(si)
95+
```
96+
97+
One could make this reponse an approximate step-function by increasing the `switch_slope` parameter, and an exact step-function by modifying expression `2` as follows.
98+
99+
```{r, eval = FALSE}
100+
alpha_t ~ alpha_max * round(sigma)
101+
```
102+
103+
However, this approach would complicate calibration by introducing a discontinuity in the log-likelihood. In practice, there is rarely a need to deviate from an approximate step function, since the approximation can be made arbitrarily sharp by increasing `switch_slope`.
104+
105+
# Simulation
106+
107+
We simulate the ODEs of this model using the following parameters.
108+
109+
```{r default}
110+
spec = mp_tmb_library("starter_models"
111+
, "si_behaviour"
112+
, package = "macpan2"
113+
)
114+
params = c("beta", "alpha_max", "threshold", "switch_slope")
115+
mp_default_list(spec)[params]
116+
```
117+
118+
119+
Simulating and plotting the state variables is done using the following code.
120+
121+
```{r simulation, fig.height=5, fig.width=4}
122+
traj = (spec
123+
|> mp_rk4()
124+
|> mp_simulator(100, c(mp_state_vars(spec), "sigma"))
125+
|> mp_trajectory()
126+
)
127+
(traj
128+
|> filter(matrix %in% mp_state_vars(spec))
129+
|> ggplot()
130+
+ aes(time, value)
131+
+ geom_line()
132+
+ facet_wrap(~matrix, scales = "free", ncol = 1)
133+
+ theme_bw()
134+
)
135+
```
136+
137+
And here are the points at which the `sigma` function was evaluated during the simulation, showing a very steep change at `I = 50` from a regime where nobody is protecting themselves, to one where all susceptible individuals protect themselves at per-capita rate `alpha_max`.
138+
139+
```{r sigmoid_plot, fig.width=4, fig.height=4}
140+
(traj
141+
|> filter(matrix %in% c("sigma", "I"))
142+
|> pivot_wider(id_cols = time, names_from = matrix)
143+
|> ggplot()
144+
+ aes(I, sigma)
145+
+ geom_point()
146+
+ theme_bw()
147+
)
148+
```
149+
150+
One can effectively make this switch a step function by increasing the `switch_slope` parameter.
151+
```{r step, fig.width=4, fig.height=4}
152+
(spec
153+
|> mp_tmb_update(default = list(switch_slope = 100))
154+
|> mp_rk4()
155+
|> mp_simulator(100, c("I", "sigma"))
156+
|> mp_trajectory()
157+
|> pivot_wider(id_cols = time, names_from = matrix)
158+
|> ggplot()
159+
+ aes(I, sigma)
160+
+ geom_point()
161+
+ theme_bw()
162+
)
163+
```

‎inst/starter_models/si_behaviour/README.md‎

Lines changed: 182 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,182 @@
1+
SI with Behaviour
2+
================
3+
Steve Walker
4+
5+
- <a href="#states" id="toc-states">States</a>
6+
- <a href="#parameters" id="toc-parameters">Parameters</a>
7+
- <a href="#dynamics" id="toc-dynamics">Dynamics</a>
8+
- <a href="#model-specification" id="toc-model-specification">Model
9+
Specification</a>
10+
- <a href="#simulation" id="toc-simulation">Simulation</a>
11+
12+
This is an SI model in which individuals can protect themselves.
13+
14+
![](./figures/flow_diagram-1.png)<!-- -->
15+
16+
The per-capita protection-rate, `alpha_t`, is time-varying, and
17+
increases as the number of infectious individuals, `I`, increases.
18+
19+
The code in this article uses the following packages.
20+
21+
``` r
22+
library(macpan2)
23+
library(ggplot2)
24+
library(dplyr)
25+
library(tidyr)
26+
```
27+
28+
This model was inspired by a question from [Irena
29+
Papst](https://github.com/papsti). The [awareness
30+
model](https://github.com/canmod/macpan2/tree/main/inst/starter_models/awareness)
31+
provides another example of a behavioural response.
32+
33+
# States
34+
35+
| Variable | Description |
36+
|----------|-----------------------------------------------------|
37+
| $S(t)$ | Number of susceptible individuals |
38+
| $I(t)$ | Number of infectious individuals |
39+
| $P(t)$ | Number of individuals who have protected themselves |
40+
41+
# Parameters
42+
43+
| Parameter | Description |
44+
|-----------------------|---------------------------------------------------------------------|
45+
| $\beta$ | Transmission rate |
46+
| $\alpha_{\text{max}}$ | Maximum per-capita protection rate |
47+
| $I^*$ | Threshold number of infectious individuals that triggers protection |
48+
| $k$ | Steepness of the sigmoidal switching function |
49+
| $N$ | Total population size (constant) |
50+
51+
# Dynamics
52+
53+
$$
54+
\begin{align*}
55+
\sigma(I) &= \frac{1}{1 + \exp\left( -k \cdot (I - I^*) \right)} \\
56+
\alpha(t) &= \alpha_{\text{max}} \cdot \sigma(I(t)) \\
57+
\\
58+
\frac{dS}{dt} &= - \beta \cdot \frac{I}{N} \cdot S - \alpha(t) \cdot S \\
59+
\frac{dI}{dt} &= \beta \cdot \frac{I}{N} \cdot S \\
60+
\frac{dP}{dt} &= \alpha(t) \cdot S
61+
\end{align*}
62+
$$
63+
64+
# Model Specification
65+
66+
This model has been specified in the `si_behaviour` directory
67+
[here](https://github.com/canmod/macpan2/blob/main/inst/starter_models/si_behaviour/tmb.R)
68+
and is accessible from the `macpan2` model library (see [Example
69+
Models](https://canmod.github.io/macpan2/articles/example_models.html)
70+
for details).
71+
72+
Printing the steps of the simulation loop illustrates how to implement a
73+
sigmoidal response of the protection rate to the prevalence of
74+
infectious individuals.
75+
76+
``` r
77+
mp_print_during(si)
78+
#> ---------------------
79+
#> At every iteration of the simulation loop (t = 1 to T):
80+
#> ---------------------
81+
#> 1: sigma ~ invlogit(switch_slope * (I - threshold))
82+
#> 2: alpha_t ~ alpha_max * sigma
83+
#> 3: mp_per_capita_flow(from = "S", to = "I", rate = "beta * I / N",
84+
#> flow_name = "infection")
85+
#> 4: mp_per_capita_flow(from = "S", to = "P", rate = "alpha_t", flow_name = "protection")
86+
```
87+
88+
One could make this reponse an approximate step-function by increasing
89+
the `switch_slope` parameter, and an exact step-function by modifying
90+
expression `2` as follows.
91+
92+
``` r
93+
alpha_t ~ alpha_max * round(sigma)
94+
```
95+
96+
However, this approach would complicate calibration by introducing a
97+
discontinuity in the log-likelihood. In practice, there is rarely a need
98+
to deviate from an approximate step function, since the approximation
99+
can be made arbitrarily sharp by increasing `switch_slope`.
100+
101+
# Simulation
102+
103+
We simulate the ODEs of this model using the following parameters.
104+
105+
``` r
106+
spec = mp_tmb_library("starter_models"
107+
, "si_behaviour"
108+
, package = "macpan2"
109+
)
110+
params = c("beta", "alpha_max", "threshold", "switch_slope")
111+
mp_default_list(spec)[params]
112+
#> $beta
113+
#> [1] 0.1
114+
#>
115+
#> $alpha_max
116+
#> [1] 0.1
117+
#>
118+
#> $threshold
119+
#> [1] 50
120+
#>
121+
#> $switch_slope
122+
#> [1] 1
123+
```
124+
125+
Simulating and plotting the state variables is done using the following
126+
code.
127+
128+
``` r
129+
traj = (spec
130+
|> mp_rk4()
131+
|> mp_simulator(100, c(mp_state_vars(spec), "sigma"))
132+
|> mp_trajectory()
133+
)
134+
(traj
135+
|> filter(matrix %in% mp_state_vars(spec))
136+
|> ggplot()
137+
+ aes(time, value)
138+
+ geom_line()
139+
+ facet_wrap(~matrix, scales = "free", ncol = 1)
140+
+ theme_bw()
141+
)
142+
```
143+
144+
![](./figures/simulation-1.png)<!-- -->
145+
146+
And here are the points at which the `sigma` function was evaluated
147+
during the simulation, showing a very steep change at `I = 50` from a
148+
regime where nobody is protecting themselves, to one where all
149+
susceptible individuals protect themselves at per-capita rate
150+
`alpha_max`.
151+
152+
``` r
153+
(traj
154+
|> filter(matrix %in% c("sigma", "I"))
155+
|> pivot_wider(id_cols = time, names_from = matrix)
156+
|> ggplot()
157+
+ aes(I, sigma)
158+
+ geom_point()
159+
+ theme_bw()
160+
)
161+
```
162+
163+
![](./figures/sigmoid_plot-1.png)<!-- -->
164+
165+
One can effectively make this switch a step function by increasing the
166+
`switch_slope` parameter.
167+
168+
``` r
169+
(spec
170+
|> mp_tmb_update(default = list(switch_slope = 100))
171+
|> mp_rk4()
172+
|> mp_simulator(100, c("I", "sigma"))
173+
|> mp_trajectory()
174+
|> pivot_wider(id_cols = time, names_from = matrix)
175+
|> ggplot()
176+
+ aes(I, sigma)
177+
+ geom_point()
178+
+ theme_bw()
179+
)
180+
```
181+
182+
![](./figures/step-1.png)<!-- -->

‎inst/starter_models/si_behaviour/figures/flow_diagram-1.png‎

10.7 KB
Loading

‎inst/starter_models/si_behaviour/figures/sigmoid_plot-1.png‎

17.4 KB
Loading

‎inst/starter_models/si_behaviour/figures/simulation-1.png‎

31.6 KB
Loading

‎inst/starter_models/si_behaviour/figures/step-1.png‎

17 KB
Loading

‎inst/starter_models/si_behaviour/tmb.R‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
library(macpan2)
2+
spec = mp_tmb_model_spec(
3+
during = list(
4+
sigma ~ invlogit(switch_slope * (I - threshold))
5+
, alpha_t ~ alpha_max * sigma
6+
, mp_per_capita_flow("S", "I", "beta * I / N", "infection")
7+
, mp_per_capita_flow("S", "P", "alpha_t", "protection")
8+
)
9+
, inits = list(S = 99, I = 1, P = 0, N = 100)
10+
, default = list(
11+
beta = 0.1 # transmission rate
12+
, alpha_max = 0.1 # maximum per-captia protection rate
13+
, threshold = 50 # threshold I, above which people protect themselves
14+
, switch_slope = 1 # larger slopes make alpha approach a step function
15+
)
16+
)

0 commit comments

Comments
 (0)