-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSIMC_SIR_Exercise_2.R
144 lines (112 loc) · 4.49 KB
/
SIMC_SIR_Exercise_2.R
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
## -----------------------------------------------------------------------------
## Please run SIMC_SIR_Exercise1.R before running this
## -----------------------------------------------------------------------------
## -----------------------------------------------------------------------------
# Initial state values
initial_state_values <- c(S = 9999,
I = 1,
R = 0
)
# Define a list of time steps
times = seq(from = 0, to = 1000, by = .1)
# Loop over a range of beta and gamma, find the size of the epidemic each time
beta.gamma.outbreak <- data.frame(beta = c(),
gamma = c(),
outbreak.size= c())
for (beta.value in seq(0,1,.1)){
for (gamma.value in seq(.1,1.,.1)){
parameters <- c(beta = beta.value, # per-contact infection rate
gamma = gamma.value # recovery rate
)
# Integrate
sir.output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters))
# Find the total outbreak size
total.recovered = sir.output[sir.output$I <1,][1,]$R
beta.gamma.outbreak = rbind(beta.gamma.outbreak,
data.frame(beta = beta.value,
gamma = gamma.value,
outbreak.size = total.recovered)
)
}
}
#beta.gamma.outbreak
## -----------------------------------------------------------------------------
p <- ggplot(data = beta.gamma.outbreak) +
geom_point(mapping = aes(x = gamma, y = beta, color = log(outbreak.size)),
size = 4) +
scale_color_gradient2() +
labs("Log(Outbreak Size)") +
scale_x_continuous(breaks = c(0,.2,.4,.6,.8,1)) +
scale_y_continuous(breaks = c(0,.2,.4,.6,.8,1)) +
xlab("gamma") +
ylab("Beta") +
theme_bw()
p
## -----------------------------------------------------------------------------
# Define a list of time steps
times = seq(from = 0, to = 1000, by = .1)
# Loop over a range of beta and gamma, find the size of the epidemic each time
r0.outbreak <- data.frame(R_0 = c(),
outbreak.size= c())
gamma.value = .2
for (beta.value in seq(0,1,.05)){
parameters <- c(beta = beta.value, # per-contact infection rate
gamma = gamma.value # recovery rate
)
# Integrate
sir.output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters))
# Find the total outbreak size
total.recovered = sir.output[sir.output$I <1,][1,]$R
r0.outbreak = rbind(r0.outbreak,
data.frame(R_0 = beta.value/gamma.value,
outbreak.size = total.recovered)
)
}
#r0.outbreak
## -----------------------------------------------------------------------------
p <- ggplot(data = r0.outbreak) +
geom_line(mapping = aes(x = R_0, y = outbreak.size), size = 2) +
xlab("R_0") +
ylab("Outbreak size") +
theme_bw()
p
## -----------------------------------------------------------------------------
# Initial state values
initial_state_values <- c(S = 9999,
I = 1,
R = 0
)
# Parameter values
parameters <- c(beta = 1, # per-contact infection rate
gamma = .2 # recovery rate
)
# Time
times = seq(from = 0, to = 50, by = .1)
# Integrate
sir.output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters))
with(as.list(parameters), {
ggplot(data = sir.output) +
geom_line(mapping = aes(x = time, y = beta/gamma * S/(S + I + R)),
size = 2 , color = 'black') +
geom_line(mapping = aes(x = time, y = I/(S+I+R)),
size = 2, color = 'red') +
xlab("Time") +
ylab("R_eff") +
theme_bw()
}
)
## -----------------------------------------------------------------------------
sir.output[sir.output$I == max(sir.output$I),]$time
## -----------------------------------------------------------------------------
with(as.list(c(sir.output[which(sir.output$I == max(sir.output$I)) ,], parameters)), {
beta/gamma * S/(S + I + R)
})