-
Notifications
You must be signed in to change notification settings - Fork 16
/
Copy pathtemporal-simulation.qmd
72 lines (61 loc) · 1.66 KB
/
temporal-simulation.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
---
title: "Simulation"
---
::: {.callout-note appearance="simple"}
This is a work in progress that is currently undergoing heavy technical editing and copy-editing
:::
```{r warning=FALSE, message=FALSE}
library(tidyverse) # essential packages
library(cowplot) # for themes
library(deSolve)
theme_set(theme_bw(base_size = 16)) # set global theme
```
## HLIR function
```{r}
HLIR_fun <- function(t, y, par) {
# Variables
H <- y[1]
L <- y[2]
I <- y[3]
R <- y[4]
beta <- par$beta
gama <- par$gama
mu <- par$mu
# Right hand side of the model
dH <- -beta * H * I
dL <- beta * H * I - gama * L
dI <- gama * L - mu * I
dR <- mu * I
return(list(c(dH, dL, dI, dR)))
}
```
## Setting up
```{r}
# Set up parameters
beta <- 0.002 # Per capita rate of infection of susceptible hosts
gama <- 1 / 5 # Rate at which exposed (i.e., latently infected) hosts become infectious
delta <- 15 # infectious period
mu <- 1 / delta
InitCond <- c(1000, 1, 0, 0)
```
## Solving equation
```{r}
steps <- seq(1, 60, by = 1)
parms <- list(beta = beta, gama = gama, mu = mu)
HLIR <- ode(InitCond, steps, HLIR_fun, parms)
epidemics <- data.frame(time = HLIR[, 1], H = HLIR[, 2], L = HLIR[, 3], I = HLIR[, 4], R = HLIR[, 5])
```
## Visualizing HLIR output
```{r}
library(ggthemes)
p1 <- epidemics %>%
ggplot() +
r4pde::theme_r4pde()+
geom_line(aes(time, H, color = "H"), size = 2) +
geom_line(aes(time, L, color = "L"), size = 2) +
geom_line(aes(time, I, color = "I"), size = 2) +
geom_line(aes(time, R, color = "R"), size = 2) +
geom_line(aes(time, I+R, color = "DIS"), size =2, linetype = 2)+
labs(y = "Population size (count)", x = "Time (days)", color = "Sites")
p1
```