-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSIMC_SIR_Exercise_4.R
127 lines (95 loc) · 4.17 KB
/
SIMC_SIR_Exercise_4.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
## -----------------------------------------------------------------------------
## Please run SIMC_SIR_Exercise1.R before running this
## -----------------------------------------------------------------------------
## -----------------------------------------------------------------------------
# Using the optimized parameter values from before
# (if you found something different, feel free to use those values here)
parameters = c(beta = 2.559449, gamma = 2.109450)
# Define times - start at week 10
times = seq(from = 10, to = 75, by = 1)
# Adjust vaccinated proportion p.vaccinated here:
p.vaccinated = .1
# Initial state values
initial_state_values <- c(S = 55000*(1 - p.vaccinated) - 1,
I = 1,
R = 55000*(p.vaccinated)
)
sir.vaccinated <- deSolve::ode(initial_state_values,
seq(from = 10, to = 75, by = 1),
sir_model,
parms = fitval$par, )
sir.vaccinated <- data.table::as.data.table(sir.vaccinated)
## -----------------------------------------------------------------------------
p <- ggplot() +
geom_line(data = sir.vaccinated,
mapping = aes(x = time, y = I),
color = 'blue', size = 2) +
geom_line(data = sir.fit,
mapping = aes(x = time, y = I),
color = 'red', size = 2 ) +
geom_point(data = plague.dat, mapping = aes(x = time, y = Deaths), size = 3) +
ylab("Deaths") + xlab("time (Weeks)") +
theme(axis.text = element_text(size = 20),
axis.title = element_text(size = 24)) +
theme_bw()
p
## -----------------------------------------------------------------------------
sir.vaccinated$R %>% tail(1)
## -----------------------------------------------------------------------------
sir.fit$R %>% tail(1)
## ----SIR model with treatment-------------------------------------------------
sir_model_treatment <- function(time, state, parameters) {
# Unpack Parameters
with(as.list(c(state,parameters)), {
# Calculate total population size
N <- S + I_t + I + R
# Calculate force of infection
lambda <- beta * (I + I_t) / N
# Calculate derivatives
dS <- -lambda * S
# Notice now that there is this new compartment of people who receive treatment
dI_t <- lambda * S * p.treatment - gamma/.75 * I_t
# And this compartment where people do not receive treatment
dI <- lambda * S * (1-p.treatment) - gamma * I
dR <- gamma * I + 2*gamma * I_t
# Return derivative
return(list(c(dS, dI, dI_t, dR)))
}
)
}
## -----------------------------------------------------------------------------
# Using the optimized parameter values from before
# (if you found something different, feel free to use those values here)
# Also specify p.treated, the fraction of infected people who receive treatment
parameters = c(beta = 2.559449, gamma = 2.109450, p.treatment = 0.25)
# Define times - start at week 10
times = seq(from = 10, to = 75, by = 1)
# Initial state values
initial_state_values <- c(S = 54999,
I = 1,
I_t = 0,
R = 0
)
sir.treated <- deSolve::ode(initial_state_values,
seq(from = 10, to = 75, by = 1),
sir_model_treatment,
parms = parameters )
sir.treated <- data.table::as.data.table(sir.treated)
## -----------------------------------------------------------------------------
p <- ggplot() +
geom_line(data = sir.treated,
mapping = aes(x = time, y = I+I_t),
color = 'blue', size = 2) +
geom_line(data = sir.fit,
mapping = aes(x = time, y = I),
color = 'red', size = 2 ) +
geom_point(data = plague.dat, mapping = aes(x = time, y = Deaths), size = 3) +
ylab("Deaths") + xlab("time (Weeks)") +
theme(axis.text = element_text(size = 20),
axis.title = element_text(size = 24)) +
theme_bw()
p
## -----------------------------------------------------------------------------
sir.treated$R %>% tail(1)
## -----------------------------------------------------------------------------
sir.fit$R %>% tail(1)