DRAFT Notes: Interventions and Calibration
Infectious Disease Modelling Specialization, Imperial College London
analysisinsightdata
Thursday, October 07, 2021
Reference: Variables, Greek Letters et al.
S = number of susceptible individuals in the population (compartment)
I = number of infected individuals in the population (compartment)
R = number of recovered/removed individuals in the population (compartment)
N = total number of individuals in the population
c = average number of contacts a susceptible makes per time
p = probability of a susceptible becoming infected on contact with infected person
beta β = infection rate (function of c * p)
gamma γ = recovery rate
lambda λ = rate of movement from exposed to infectious (beta * I/N)
R0 R0 R-naught = basic reproduction number of infectious agents (beta / gamma)
Reff Reff Rt = Effective Reproduction Number (R0 * S/N)
pvacc p = vaccination rate or coverage
pvacc_threshold pc = critical vaccination threshold; (1 - (1/R0)); or (1 - (1/R0))/(1 - cs); or (1 - (1/R0))/(1 - (cs * ci))
veff veff = vaccine efficacy; reduces susceptibility by (1 - cs); reduces infection by (1 - ci)
peff peff = effective vaccine coverage (in all-or-nothing = vaccination rate p * veff vaccine efficacy)
cs c_s cS = effectiveness of the vaccine in reducing susceptibility of vaccinated individuals; vaccination reduces susceptibility by (1 - cs)
ci c_i ci = effectiveness of the vaccine in reducing infectivity of vaccinated individuals; vaccination reduces infectivity by (1 - ci)
sigma σ = waning immunity rate (some references uses xi ξ)
mu μ = background mortality rate or infection-related excess mortality rate (depending on model)
birth b = immigration or birth rate (same as mu if stable population)
M = number of deceased (compartment)
T = separate compartment for infected people on treatment; if used, I becomes the compartment for untreated infected people
h = rate of transition from untreated (I) infected people to treated (T) infected people
gammaT γT = recovery rate from treated (T) infected people
betaT βT = infection rate for treated (T) infected people
V = number of vaccinated people; can be used to track vaccine-induced immunity and model leaky scenarios (compartment)
Iv Iv = number of infected vaccinated individuals (compartment)
I0 I0 I-naught = initial prevalence at time t equals 0
E = separate compartment for exposed people in the latent period between being infected and becoming infectious (infected but not yet infectious)
delta δ = rate of progression to infectious from exposed (coursework uses alpha α, other references use sigma σ)
Modelling treatment
Extend the simple SIR model to incorporate a treatment intervention and study its effect on the epidemiology of an immunising infection.
Model a viral epidemic in a town of 300,000 people, which starts with the introduction of a single infected case in a fully susceptible population. Public health authorities believe the infection rate is 0.6 days\(^{-1}\) and that people stay infected for 5 days on average before recovering. Once recovered, people are immune to reinfection.
A treatment against the virus is available. Although it does not cure people immediately, it helps them to recover more quickly: on average, people recover after 1.25 days of treatment. Treatment has no effect on transmission of the virus, so treated people are just as infectious as untreated infected people. Imagine that people start taking the treatment on average 4 days after becoming infected.
Extend the simple SIR model to include a separate compartment for people on treatment, \(T\). Here, we have a new situation where infectious people are now separated into 2 different compartments: treated (\(T\)) and untreated (\(I\)) infected people can transmit the virus. This means that we have to adjust the force of infection \(\lambda\), as described below.
Specifying the force of infection in the treatment model:
The force of infection \(\lambda\) acts on susceptibles. It removes people from the \(S\) compartment and adds them into the \(I\) compartment - this remains the same in our treatment model.
The only difference is that people in both the \(I\) and \(T\) compartments are a source of infection now. Here, we assume they are equally good at passing on the infection, which means they both transmit the virus to susceptibles at a rate \(\beta\). \(\beta\) is simply the average number of secondary infections caused by an infectious person per unit time - here we are assuming that untreated and treated people cause the same average number of secondary infections per unit time.
In the basic SIR model, the equation for the force of infection was:
\[\begin{align} \lambda & = \beta \frac{I}{N} \end{align}\]
nicesubtitle <- "SIR Model v1 basic model"
print(nicesubtitle)## [1] "SIR Model v1 basic model"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
duration <- 365 # total number of days
tsteps <- 1 # chunk into 1 day intervals
beta <- 0.6 # infection rate given as 1 day^-1
gamma <- 1/5 # recover after 5 days for untreated I
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate for untreated I
R0 = R0
)) ## beta gamma R0
## 0.6 0.2 3.0initial_state_values <- c(S = N-1,
I = 1,
R = 0)
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dR <- (gamma * I)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R),digits=5)*100,
preval_Inf = round(I/(S+ I + R),digits=5)*100,
propor_Re = round(R/(S + I + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output %>%
arrange(-I, time) %>%
head(1)## time S I R still_Su preval_Inf propor_Re Reff
## 1 33 108507.7 89796.58 101695.7 36.169 29.932 33.899 1.085077print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S I R still_Su preval_Inf propor_Re Reff
## 1 34 90632.92 89671.04 119696 30.211 29.89 39.899 0.9063292print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I R still_Su preval_Inf propor_Re Reff
## 366 365 17855.96 2.732142e-18 282144 5.952 0 94.048 0.1785596print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output %>%
select(time, Reff) %>%
filter(time < 120) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
All we have to do now is add those additional infectious people in the \(T\) compartment:
\[\begin{align} \lambda & = \beta \frac{I+T}{N} \end{align}\]
Or another way of expressing this:
\[\begin{align} \lambda & = \beta \frac{I}{N} + \beta \frac{T}{N} \end{align}\]
Investigate the impact that treatment has on the course of this epidemic. Use the information above to define a new model function, initial conditions and parameters and solve this new system.
Plot the prevalence of infection over time, being careful to include all infected compartments.
Question: How many people are infected at the peak of the epidemic?
- There were 37,615 (12.5%) infected at the peak on day 42.
nicesubtitle <- "SIR Model v2a with untreated I and treated T infected"
print(nicesubtitle)## [1] "SIR Model v2a with untreated I and treated T infected"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
duration <- 365 # total number of days
tsteps <- 1 # chunk into 1 day intervals
beta <- 0.6 # infection rate given as 1 day^-1
gamma <- 1/5 # recover after 5 days for untreated I
gammaT <- 1/1.25 # recover after 1.25 days for treated T
h <- 1/4 # rate of transition from I to T
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate for untreated I
gammaT = gammaT, # recovery rate for treated T
h = h, # rate of transition from I to T
R0 = R0
)) ## beta gamma gammaT h R0
## 0.60 0.20 0.80 0.25 3.00initial_state_values <- c(S = N-1,
I = 1,
T = 0,
R = 0)
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_c2a <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + T + R
lambda <- beta * (I + T)/N
dS <- -(lambda * S)
dI <- (lambda * S) -((gamma + h) * I)
dT <- (h * I) -(gammaT * T)
dR <- (gamma * I) + (gammaT * T)
return(list(c(dS, dI, dT, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_c2a,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + T + R),digits=5)*100,
preval_Inf = round(I/(S+ I + T + R),digits=5)*100,
in_Treated = round(T/(S+ I + T + R),digits=5)*100,
propor_Re = round(R/(S + I + T + R),digits=5)*100,
combo_I_T = round((I+T),digits=1),
perc_I_T = round((I+T)/(S + I + T + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + T + R)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output %>%
arrange(-I, time) %>%
head(1)## time S I T R still_Su preval_Inf in_Treated
## 1 42 169285.2 28881.29 8733.865 93099.67 56.428 9.627 2.911
## propor_Re combo_I_T perc_I_T Reff
## 1 31.033 37615.2 12.538 1.692852print("peak infection & treated days when (I + T) is at its max: ")## [1] "peak infection & treated days when (I + T) is at its max: "output %>%
arrange(-(I+T), time) %>%
head(1)## time S I T R still_Su preval_Inf in_Treated
## 1 42 169285.2 28881.29 8733.865 93099.67 56.428 9.627 2.911
## propor_Re combo_I_T perc_I_T Reff
## 1 31.033 37615.2 12.538 1.692852print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S I T R still_Su preval_Inf in_Treated
## 1 52 99060.75 9410.079 3683.139 187846 33.02 3.137 1.228
## propor_Re combo_I_T perc_I_T Reff
## 1 62.615 13093.2 4.364 0.9906075print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I T R still_Su preval_Inf
## 366 365 86189.98 4.096914e-23 1.721887e-23 213810 28.73 0
## in_Treated propor_Re combo_I_T perc_I_T Reff
## 366 0 71.27 0 0 0.8618998print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -in_Treated, -propor_Re, -Reff,
-combo_I_T, -perc_I_T) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","purple","green")) +
scale_shape_manual(values = c(0,4,7,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" h of", parameters['h'],
" gammaT of", parameters['gammaT'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output %>%
select(time, Reff) %>%
filter(time < 120) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" h of", parameters['h'],
" gammaT of", parameters['gammaT'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
Despite the treatment, an epidemic still occurs, albeit a smaller one than if no one was treated.
How you could improve the impact of this treatment on reducing the prevalence of this infection, by changing the treatment initiation rate?
- Changing the transition rate from 1/4 to 1 (treat everyone after 1 day) reduces peak infected (I + T) to 2,362 (0.8%), pushing out to day 111. Although we have no infected after 365 days, We still need to bump up the transition rate to treatment (h) such that there is no epidemic.
nicesubtitle <- "SIR Model v2a2 with untreated I and treated T infected"
print(nicesubtitle)## [1] "SIR Model v2a2 with untreated I and treated T infected"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
duration <- 365 # total number of days
tsteps <- 1 # chunk into 1 day intervals
beta <- 0.6 # infection rate given as 1 day^-1
gamma <- 1/5 # recover after 5 days for untreated I
gammaT <- 1/1.25 # recover after 1.25 days for treated T
h <- 1 # rate of transition from I to T
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate for untreated I
gammaT = gammaT, # recovery rate for treated T
h = h, # rate of transition from I to T
R0 = R0
)) ## beta gamma gammaT h R0
## 0.6 0.2 0.8 1.0 3.0initial_state_values <- c(S = N-1,
I = 1,
T = 0,
R = 0)
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_c2a <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + T + R
lambda <- beta * (I + T)/N
dS <- -(lambda * S)
dI <- (lambda * S) -((gamma + h) * I)
dT <- (h * I) -(gammaT * T)
dR <- (gamma * I) + (gammaT * T)
return(list(c(dS, dI, dT, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_c2a,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + T + R),digits=5)*100,
preval_Inf = round(I/(S+ I + T + R),digits=5)*100,
in_Treated = round(T/(S+ I + T + R),digits=5)*100,
propor_Re = round(R/(S + I + T + R),digits=5)*100,
combo_I_T = round((I+T),digits=1),
perc_I_T = round((I+T)/(S + I + T + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + T + R)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output %>%
arrange(-I, time) %>%
head(1)## time S I T R still_Su preval_Inf in_Treated
## 1 110 267997.9 1051.608 1305.755 29644.78 89.333 0.351 0.435
## propor_Re combo_I_T perc_I_T Reff
## 1 9.882 2357.4 0.786 2.679979print("peak infection & treated days when (I + T) is at its max: ")## [1] "peak infection & treated days when (I + T) is at its max: "output %>%
arrange(-(I+T), time) %>%
head(1)## time S I T R still_Su preval_Inf in_Treated
## 1 111 266735.6 1051.58 1310.75 30902.02 88.912 0.351 0.437
## propor_Re combo_I_T perc_I_T Reff
## 1 10.301 2362.3 0.787 2.667356print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## [1] time S I T R still_Su
## [7] preval_Inf in_Treated propor_Re combo_I_T perc_I_T Reff
## <0 rows> (or 0-length row.names)print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I T R still_Su preval_Inf
## 366 365 235892.1 1.477466e-05 2.04156e-05 64107.89 78.631 0
## in_Treated propor_Re combo_I_T perc_I_T Reff
## 366 0 21.369 0 0 2.358921print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -in_Treated, -propor_Re, -Reff,
-combo_I_T, -perc_I_T) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","purple","green")) +
scale_shape_manual(values = c(0,4,7,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" h of", parameters['h'],
" gammaT of", parameters['gammaT'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output %>%
select(time, Reff) %>%
filter(time < 120) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" h of", parameters['h'],
" gammaT of", parameters['gammaT'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
Question: How rapidly does treatment need to be initiated in order to interrupt transmission, i.e. to bring R0 below 1? Based on this, do you think it is feasible to interrupt transmission through treatment alone?
To interrupt transmission by bringing R0 below 1, the treatment initiation rate needs to be at least 1.6 per day, which means people need to start treatment less than a day after becoming infected on average.
“To achieve a reduction in the treatment initiation rate, think about what it depends on: the time it takes people to go to a doctor, the time it takes to get a diagnosis, and the time from diagnosis to starting treatment for example. These in turn depend on many situational aspects such as whether the disease is symptomatic, the healthcare system, which test is required for a diagnosis etc. One way of increasing the rate of treatment initiation would be through active case-finding for example, rather than waiting for people to seek medical attention themselves.”
“However, given all these factors, achieving a treatment initiation rate as high as required in this example does not seem feasible.”
It is also possible to derive a formula to calculate R0 for more complex extensions to the SIR model like this one. Derive and apply this equation to your parameter values to determine the minimum treatment initiation rate that you would need to reduce R0 below 1. For this, remember the definition of R0.
Question: Using only reasoning based on R0, what is the minimum value of h needed to interrupt transmission?
R0 is defined as the average number of secondary infections caused by a single infectious case (in a totally susceptible population). We can derive the following equation:
\[\begin{align} R_0 & = \frac{\beta}{\gamma+h} + \frac{h}{\gamma+h} \times \frac{\beta}{\gamma_T} \end{align}\]
Here, we are taking the average of secondary infections caused by index cases in the I and the T compartment, keeping in mind that only a proportion \[\begin{align} \frac{h}{\gamma+h}\end{align}\] move into the treatment compartment before recovering.
Solving this equation with our parameter values to obtain R0 < 1:
\[\begin{align} 1 & > \frac{0.6}{0.2+h} + \frac{h}{0.2+h} \times \frac{0.6}{0.8} \\ \\ 0.2 + h & > 0.6 + 0.75 h \\ \\ h & > 1.6 \end{align}\]
- Changing the transition rate to 1.6 (treat everyone after about half a day – which is unrealistic) reduces peak infected (I + T) to 1.2 (about 0%), on day 5. Since we started with 1 infected person, this pretty much equates to no epidemic.
nicesubtitle <- "SIR Model v2a3 with untreated I and treated T infected"
print(nicesubtitle)## [1] "SIR Model v2a3 with untreated I and treated T infected"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
duration <- 365 # total number of days
tsteps <- 1 # chunk into 1 day intervals
beta <- 0.6 # infection rate given as 1 day^-1
gamma <- 1/5 # recover after 5 days for untreated I
gammaT <- 1/1.25 # recover after 1.25 days for treated T
h <- 1.6 # rate of transition from I to T
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate for untreated I
gammaT = gammaT, # recovery rate for treated T
h = h, # rate of transition from I to T
R0 = R0
)) ## beta gamma gammaT h R0
## 0.6 0.2 0.8 1.6 3.0initial_state_values <- c(S = N-1,
I = 1,
T = 0,
R = 0)
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_c2a <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + T + R
lambda <- beta * (I + T)/N
dS <- -(lambda * S)
dI <- (lambda * S) -((gamma + h) * I)
dT <- (h * I) -(gammaT * T)
dR <- (gamma * I) + (gammaT * T)
return(list(c(dS, dI, dT, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_c2a,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + T + R),digits=5)*100,
preval_Inf = round(I/(S+ I + T + R),digits=5)*100,
in_Treated = round(T/(S+ I + T + R),digits=5)*100,
propor_Re = round(R/(S + I + T + R),digits=5)*100,
combo_I_T = round((I+T),digits=1),
perc_I_T = round((I+T)/(S + I + T + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + T + R)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output %>%
arrange(-I, time) %>%
head(1)## time S I T R still_Su preval_Inf in_Treated propor_Re combo_I_T perc_I_T
## 1 0 299999 1 0 0 100 0 0 0 1 0
## Reff
## 1 2.99999print("peak infection & treated days when (I + T) is at its max: ")## [1] "peak infection & treated days when (I + T) is at its max: "output %>%
arrange(-(I+T), time) %>%
head(1)## time S I T R still_Su preval_Inf in_Treated
## 1 5 299995.5 0.4000112 0.7999415 3.339973 99.998 0 0
## propor_Re combo_I_T perc_I_T Reff
## 1 0.001 1.2 0 2.999955print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## [1] time S I T R still_Su
## [7] preval_Inf in_Treated propor_Re combo_I_T perc_I_T Reff
## <0 rows> (or 0-length row.names)print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I T R still_Su preval_Inf in_Treated
## 366 365 299746.1 0.3568364 0.7142162 252.8295 99.915 0 0
## propor_Re combo_I_T perc_I_T Reff
## 366 0.084 1.1 0 2.997461print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -in_Treated, -propor_Re, -Reff,
-combo_I_T, -perc_I_T) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","purple","green")) +
scale_shape_manual(values = c(0,4,7,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" h of", parameters['h'],
" gammaT of", parameters['gammaT'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output %>%
select(time, Reff) %>%
filter(time < 120) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" h of", parameters['h'],
" gammaT of", parameters['gammaT'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
Question: What other (theoretical) changes could you make to this treatment to improve its impact on the epidemic?
- On top of quickly responding with treatment, perhaps using a stronger treatment dose or an adjuvant might possibly help increase the recovery rate [from treatment], or decrease its ability to spread. Public health measures (isolation / quarantine, for example) also supports this progress.
A separate compartment for vaccination
Develop a simple model of vaccination that you can modify to explore different ways of modelling vaccine impact. Previously, we represented vaccination by choosing a certain proportion of population immunity, represented by the \(R\) compartment, at the beginning of the simulation. If we want to distinguish the effect of vaccination, or if we want to capture some important difference in vaccine-induced immunity compared to recovery-induced immunity, we can include a separate compartment \(V\) for vaccinated people in addition to the recovered compartment \(R\). In this exercise we are still modelling vaccination coverage through the initial conditions, so no new parameters are added to the model, but the vaccination coverage p should be included in the initial state values. The proportion of the population that is effectively vaccinated does not change over the simulation.
Draw a diagram of the model structure and write down the initial conditions for a scenario of a single infected person introduced into an unexposed population with pre-existing vaccination at coverage p.
\[\begin{align} S0 & = (1-p) N \\ I0 & = 1 \\ R0 & = 0 \\ V0 & = p N \end{align}\]
where p [coded here as pvacc] is the effective vaccination coverage and \(N\) is the total population size.
Reuse the code from your simple SIR model and add this new compartment in the cell below. Use model parameters so that R0 > 1 and total population size assumptions of your choice to double-check that your code works, by plotting the proportion of the population in each compartment over time.
Three easy-to-check conditions that your output needs to fulfill:
- With no or a low vaccination coverage, you should see an epidemic (given that R0 > 1).
- With a R0 of 3.0 and a pvacc of 0.10, infection peaks at Day 38 with 70,463 (23.5%) people.
- With a vaccination coverage > herd immunity threshold, you should see that no epidemic occurs.
- With a R0 of 3.0, we calculate pc pvacc_thresh of 0.67. Using that value as the initial vaccination coverage p pvacc, we observe no epidemic.
- The proportion of the population in the vaccinated compartment V should stay the same over time.
- Vaccinated numbers in this model remains static.
nicesubtitle <- "SIR Model v3a with separate V vaccinated"
print(nicesubtitle)## [1] "SIR Model v3a with separate V vaccinated"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
duration <- 120 # total number of days
tsteps <- 1 # chunk into 1 day intervals
beta <- 0.6 # infection rate given as 1 day^-1
gamma <- 1/5 # recover after 5 days for untreated I
pvacc <- .10 # vaccination coverage aka p
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0)
(parameters <- c(
beta = beta,
gamma = gamma,
R0 = R0,
pvacc = pvacc,
pvacc_thresh = pvacc_thresh
))## beta gamma R0 pvacc pvacc_thresh
## 0.6000000 0.2000000 3.0000000 0.1000000 0.6666667initial_state_values <- c(S = round((N-1) * (1-pvacc)),
I = 1,
R = 0,
V = round((N-1) * (pvacc)) )
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_c3a <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R + V
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dR <- (gamma * I)
dV <- 0 # V stays the same over time,
# so rate of change equals 0
return(list(c(dS, dI, dR, dV)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_c3a,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R + V),digits=5)*100,
preval_Inf = round(I/(S+ I + R + V),digits=5)*100,
propor_Re = round(R/(S + I + R + V),digits=5)*100,
propor_Vac = round(V/(S + I + R + V),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R + V)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output %>%
arrange(-I, time) %>%
head(1)## time S I R V still_Su preval_Inf propor_Re
## 1 38 106655.7 70463.08 92881.21 30000 35.552 23.488 30.96
## propor_Vac Reff
## 1 10 1.066557print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S I R V still_Su preval_Inf propor_Re
## 1 39 92612.54 70388.12 106999.3 30000 30.871 23.463 35.666
## propor_Vac Reff
## 1 10 0.9261254print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I R V still_Su preval_Inf propor_Re
## 121 120 22790.05 0.5546693 247209.4 30000 7.597 0 82.403
## propor_Vac Reff
## 121 10 0.2279005print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -propor_Vac, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green","steelblue")) +
scale_shape_manual(values = c(0,4,1,20)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output %>%
select(time, Reff) %>%
filter(time < 120) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
- Setting initial vaccination coverage p pvacc equal to calculate pc pvacc_thresh of 0.67.
nicesubtitle <- "SIR Model v3a2 with separate V vaccinated"
print(nicesubtitle)## [1] "SIR Model v3a2 with separate V vaccinated"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
duration <- 120 # total number of days
tsteps <- 1 # chunk into 1 day intervals
beta <- 0.6 # infection rate given as 1 day^-1
gamma <- 1/5 # recover after 5 days for untreated I
pvacc <- .10 # vaccination coverage aka p
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0)
pvacc <- pvacc_thresh # set vaccination coverage to cvt
(parameters <- c(
beta = beta,
gamma = gamma,
R0 = R0,
pvacc = pvacc,
pvacc_thresh = pvacc_thresh
))## beta gamma R0 pvacc pvacc_thresh
## 0.6000000 0.2000000 3.0000000 0.6666667 0.6666667initial_state_values <- c(S = round((N-1) * (1-pvacc)),
I = 1,
R = 0,
V = round((N-1) * (pvacc)) )
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_c3a <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R + V
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dR <- (gamma * I)
dV <- 0 # V stays the same over time,
# so rate of change equals 0
return(list(c(dS, dI, dR, dV)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_c3a,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R + V),digits=5)*100,
preval_Inf = round(I/(S+ I + R + V),digits=5)*100,
propor_Re = round(R/(S + I + R + V),digits=5)*100,
propor_Vac = round(V/(S + I + R + V),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R + V)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output %>%
arrange(-I, time) %>%
head(1)## time S I R V still_Su preval_Inf propor_Re propor_Vac Reff
## 1 0 1e+05 1 0 199999 33.333 0 0 66.666 1print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S I R V still_Su preval_Inf propor_Re propor_Vac Reff
## 1 0 1e+05 1 0 199999 33.333 0 0 66.666 1print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I R V still_Su preval_Inf propor_Re
## 121 120 99976.03 0.9971257 23.97699 199999 33.325 0 0.008
## propor_Vac Reff
## 121 66.666 0.9997603print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -propor_Vac, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green","steelblue")) +
scale_shape_manual(values = c(0,4,1,20)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output %>%
select(time, Reff) %>%
filter(time < 120) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
Modelling a leaky vaccine
We will continue working on the model that assumed that a proportion p of the population received a vaccine that is fully effective in protecting against infection.
Question: Assuming \(\beta\) equals 0.25 days\(^{-1}\) and \(\gamma\) equals 0.1 days\(^{-1}\), what proportion of the population would have to be vaccinated with a perfectly effective vaccine to prevent an epidemic?
- It appears that pvacc p should be set to the pvacc_thresh pc critical vaccination threshold of 0.60 to avoid an epidemic.
Using the formula for the herd immunity threshold, we need a vaccine coverage of 60% with a perfect vaccine:
\[\begin{align} p_c & = 1- \frac{1}{R_{0}} \\ & = 1- \frac{\gamma}{\beta} \\ & = 1 - \frac{0.1}{0.25} \\ & = 0.6 \end{align}\]
nicesubtitle <- "SIR Model v3a3 with perfectly eff vaccine"
print(nicesubtitle)## [1] "SIR Model v3a3 with perfectly eff vaccine"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
duration <- 720 # total number of days
tsteps <- 1 # chunk into 1 day intervals
beta <- 0.25 # infection rate given as 1 day^-1
gamma <- 0.1 # recover after 5 days for untreated I
pvacc <- .10 # vaccination coverage aka p
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0)
pvacc <- pvacc_thresh # set vaccination coverage to cvt
(parameters <- c(
beta = beta,
gamma = gamma,
R0 = R0,
pvacc = pvacc,
pvacc_thresh = pvacc_thresh
))## beta gamma R0 pvacc pvacc_thresh
## 0.25 0.10 2.50 0.60 0.60initial_state_values <- c(S = round((N-1) * (1-pvacc)),
I = 1,
R = 0,
V = round((N-1) * (pvacc)) )
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_c3a <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R + V
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dR <- (gamma * I)
dV <- 0 # V stays the same over time,
# so rate of change equals 0
return(list(c(dS, dI, dR, dV)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_c3a,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R + V),digits=5)*100,
preval_Inf = round(I/(S+ I + R + V),digits=5)*100,
propor_Re = round(R/(S + I + R + V),digits=5)*100,
propor_Vac = round(V/(S + I + R + V),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R + V)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output %>%
arrange(-I, time) %>%
head(1)## time S I R V still_Su preval_Inf propor_Re propor_Vac Reff
## 1 0 120000 1 0 179999 40 0 0 60 1print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S I R V still_Su preval_Inf propor_Re propor_Vac Reff
## 1 1 119999.9 1 0.1 179999 40 0 0 60 0.9999992print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I R V still_Su preval_Inf propor_Re
## 721 720 119928.5 0.9787115 71.48612 179999 39.976 0 0.024
## propor_Vac Reff
## 721 60 0.9994045print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -propor_Vac, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green","steelblue")) +
scale_shape_manual(values = c(0,4,1,20)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output %>%
select(time, Reff) %>%
filter(time < 120) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
What if the vaccine is not perfect (as is usually the case in reality)? There are two common assumptions modellers make about how imperfect vaccines work: the “all-or-nothing” and the “leaky” mechanism. All-or-nothing vaccines give full immunity to a proportion of those vaccinated, but don’t work at all in the other vaccinated people.
Question: Given the parameter assumptions above, what proportion of the population would have to be vaccinated with an all-or-nothing vaccine with 70% efficacy to prevent an epidemic?
“Under the assumption of an all-or-nothing vaccine, we can simply multiply the vaccine efficacy veff and the vaccine coverage to calculate the effective coverage peff:
\[\begin{align} v_{eff} * p_{eff} & = 0.6 \\ p_{eff} = \frac{0.6}{0.7} \\ p_{eff} = 0.86 \end{align}\]
Therefore, we need at least 86% coverage of a leaky vaccine with efficacy of 70%, to interrupt transmission (R0 < 1).”
In contrast, leaky vaccines reduce susceptibility in everyone who received the vaccine by the same fraction. This means that vaccinated people can still become infected, but at a reduced rate (force of infection), like so:
Question: Based on the diagram, write out the differential equations for this model. What is the value of cS for a leaky vaccine with 70% efficacy?
” \[\begin{align} \frac{dS}{dt} & = -\beta \frac{I}{N} S \\ \frac{dI}{dt} & = \beta \frac{I}{N} S + c_{s} \beta \frac{I}{N} V - \gamma I \\ \frac{dR}{dt} & = \gamma I \\ \frac{dV}{dt} & = -c_{s} \beta \frac{I}{N} V \end{align}\]
For a leaky vaccine with 70% efficacy, the value of cs would be 0.3, reflecting the degree to which susceptibility is reduced.”
Building on where you set up the \(V\) compartment for vaccinated people, explore the impact of a “leaky” vaccine on an epidemic. Adapt the differential equations to reflect the leaky mechanism. Vaccination still happens before before the epidemic, so represent vaccination coverage through the initial conditions.
Question: With a leaky vaccine with 70% efficacy, what proportion of the population would have to be vaccinated to prevent an epidemic (Reff < 1)?
Using vaccination coverage p = .60 for a similar vaccine but with only 70% efficacy does not prevent an epidemic. As seen in the calculations above, we need to factor the 70% vaccination efficacy into the critical vaccination threshold calculation, which results in a cvt of .86 to prevent an epidemic.
- Using vaccination coverage p = .60 for a vaccine with 70% efficacy does not prevent an epidemic:
nicesubtitle <- "SIR Model v4a1 leaky vaccine w/ 70% efficacy"
print(nicesubtitle)## [1] "SIR Model v4a1 leaky vaccine w/ 70% efficacy"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
duration <- 720 # total number of days
tsteps <- 1 # chunk into 1 day intervals
veff <- .70 # vaccine efficacy
cs <- 1 - veff # reduction in the force of infection (1 - % eff)
beta <- 0.25 # infection rate given as 1 day^-1
gamma <- 0.1 # recovery rate given as 1 day^-1
pvacc <- .60 # vaccination coverage aka p
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0) # threshold for vaccine 100% efficacy
# pvacc <- pvacc_thresh # set vaccination coverage to cvt
peff <- pvacc * veff # effective vaccine coverage
pvacc_thresh <- (1 - (1/R0))/veff # recalc results to consider efficacy
# pvacc <- pvacc_thresh # set vaccination coverage to recalc cvt
(parameters <- c(
veff = veff,
peff = peff,
cs = cs,
beta = beta,
gamma = gamma,
R0 = R0,
pvacc = pvacc,
pvacc_thresh = pvacc_thresh
))## veff peff cs beta gamma R0
## 0.7000000 0.4200000 0.3000000 0.2500000 0.1000000 2.5000000
## pvacc pvacc_thresh
## 0.6000000 0.8571429initial_state_values <- c(S = round((N-1) * (1-pvacc)),
I = 1,
R = 0,
V = round((N-1) * (pvacc)) )
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_c4a <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R + V
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I) +(cs * (lambda * V))
dR <- (gamma * I)
dV <- -(cs * (lambda * V)) # leaky vacc only 70% eff
return(list(c(dS, dI, dR, dV)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_c4a,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R + V),digits=5)*100,
preval_Inf = round(I/(S+ I + R + V),digits=5)*100,
propor_Re = round(R/(S + I + R + V),digits=5)*100,
propor_Vac = round(V/(S + I + R + V),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R + V)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output %>%
arrange(-I, time) %>%
head(1)## time S I R V still_Su preval_Inf propor_Re
## 1 245 73334.63 12292.84 59095.06 155277.5 24.445 4.098 19.698
## propor_Vac Reff
## 1 51.759 0.6111219print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S I R V still_Su preval_Inf propor_Re
## 1 1 119999.9 1.046028 0.1022851 179999 40 0 0
## propor_Vac Reff
## 1 60 0.9999991print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I R V still_Su preval_Inf propor_Re
## 721 720 41303.23 0.006452345 127986.1 130710.6 13.768 0 42.662
## propor_Vac Reff
## 721 43.57 0.3441935print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -propor_Vac, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green","steelblue")) +
scale_shape_manual(values = c(0,4,1,20)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" cs of", parameters['cs'],
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output %>%
select(time, Reff) %>%
filter(time < 120) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" cs of", parameters['cs'],
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
- Using critical vaccination threshold calculated with vaccination efficacy of 70% prevents an epidemic:
nicesubtitle <- "SIR Model v4a2 leaky vaccine w/ 70% efficacy"
print(nicesubtitle)## [1] "SIR Model v4a2 leaky vaccine w/ 70% efficacy"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 300000 # population size
duration <- 720 # total number of days
tsteps <- 1 # chunk into 1 day intervals
veff <- .70 # vaccine efficacy
cs <- 1 - veff # reduction in the force of infection (1 - % eff)
beta <- 0.25 # infection rate given as 1 day^-1
gamma <- 0.1 # recovery rate given as 1 day^-1
pvacc <- .60 # vaccination coverage aka p
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0) # threshold for vaccine 100% efficacy
# pvacc <- pvacc_thresh # set vaccination coverage to cvt
peff <- pvacc * veff # effective vaccine coverage
pvacc_thresh <- (1 - (1/R0))/veff # recalc results to consider efficacy
pvacc <- pvacc_thresh # set vaccination coverage to recalc cvt
(parameters <- c(
veff = veff,
peff = peff,
cs = cs,
beta = beta,
gamma = gamma,
R0 = R0,
pvacc = pvacc,
pvacc_thresh = pvacc_thresh
))## veff peff cs beta gamma R0
## 0.7000000 0.4200000 0.3000000 0.2500000 0.1000000 2.5000000
## pvacc pvacc_thresh
## 0.8571429 0.8571429initial_state_values <- c(S = round((N-1) * (1-pvacc)),
I = 1,
R = 0,
V = round((N-1) * (pvacc)) )
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_c4a <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R + V
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I) +(cs * (lambda * V))
dR <- (gamma * I)
dV <- -(cs * (lambda * V)) # leaky vacc only 70% eff
return(list(c(dS, dI, dR, dV)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_c4a,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R + V),digits=5)*100,
preval_Inf = round(I/(S+ I + R + V),digits=5)*100,
propor_Re = round(R/(S + I + R + V),digits=5)*100,
propor_Vac = round(V/(S + I + R + V),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R + V)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output %>%
arrange(-I, time) %>%
head(1)## time S I R V still_Su preval_Inf propor_Re propor_Vac Reff
## 1 0 42857 1 0 257142 14.286 0 0 85.714 0.3571417print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S I R V still_Su preval_Inf propor_Re propor_Vac Reff
## 1 0 42857 1 0 257142 14.286 0 0 85.714 0.3571417print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I R V still_Su preval_Inf propor_Re
## 721 720 42831.4 0.9879791 71.70769 257095.9 14.277 0 0.024
## propor_Vac Reff
## 721 85.699 0.3569283print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -propor_Re, -propor_Vac, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green","steelblue")) +
scale_shape_manual(values = c(0,4,1,20)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" cs of", parameters['cs'],
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output %>%
select(time, Reff) %>%
filter(time < 120) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" cs of", parameters['cs'],
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
Question: Confirm this result using equations.
“We can calculate the critical vaccination coverage needed to interrupt transmission (Reff < 1), in the following way.
Remember that in a simple homogenous model, Reff is proportional to the number of susceptible people in the population. In this case:
\[\begin{align} R_{eff} & = (1-p) \times R_0 + p c_S \times R_0 \end{align}\]
where p is the proportion of the population receiving the vaccine, and cs is the reduction in susceptibility owing to the vaccine.
Setting p = pc when Reff = 1, and solving this to find pc gives:
\[\begin{align} p_c & = \frac{1-\frac{1}{R_0}}{1-c_s} \\ p_c & = \frac{1-\frac{1}{\frac{0.25}{0.1}}}{1-0.3} \\ p_c & = 0.86 \end{align}\] ”
Modelling additional vaccine effects
Model a vaccine with a combined leaky effect. The vaccine reduces not only susceptibility, but also the infectivity of those who become infected despite being vaccinated. The parameter cS describes the effectiveness of the vaccine in reducing susceptibility of vaccinated individuals, and the parameter ci describes the effectiveness of the vaccine in reducing infectivity of vaccinated individuals.
Differential equations for a model that describes the impact of a vaccine with a combined leaky mechanism.
”
\[\begin{align}
\frac{dS}{dt} & = -(\beta \frac{I}{N}+ c_{i} \beta \frac{I_{V}}{N}) S \\
\\
\frac{dI}{dt} & = (\beta \frac{I}{N}+ c_{i} \beta \frac{I_{V}}{N}) S - \gamma I \\
\\
\frac{dI_{V}}{dt} & = c_{s} (\beta \frac{I}{N} + c_{i} \beta \frac{I_{V}}{N}) V - \gamma I_{V} \\
\\
\frac{dR}{dt} & = \gamma I + \gamma I_{V} \\
\\
\frac{dV}{dt} & = - c_{s}(\beta \frac{I}{N} + c_{i} \beta \frac{I_{V}}{N}) V
\end{align}\] ”
Incorporate the new compartment and parameter with the same parameter values as in the previous activity. You can assume a single infected person is introduced into an otherwise susceptible population. Code the force of infection for a model with more than one infectious compartment.
Question: What is the efficacy of the vaccine in terms of reducing susceptibility and in terms of reducing infectivity?
- A cs value of .3 equates to a (1 - cs =) 70% effective vaccination, reducing susceptibility by 70%
- A ci value of .5 means that this vaccination reduces infectivity by (1 - ci =) 50%
“Vaccine efficacy in terms of reducing susceptibility is \((1-c_s) \times 100\) = 70% and the vaccine efficacy in terms of reducing infectivity is \((1-c_i) \times 100\) = 50%. This means that those who are infected after being vaccinated are half as infectious as those who never received the vaccine.”
Question: What is the peak prevalence (peak number of infected people) with a vaccine coverage of 60%?
- At day 478, the value of (I + Iv) is 12,304 or 1.23% of the population
- Using vaccination coverage p = .60 does not prevent an epidemic.
nicesubtitle <- "SIR Model v4b1 vaccine w/ combined leaky effects"
print(nicesubtitle)## [1] "SIR Model v4b1 vaccine w/ combined leaky effects"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 1000000 # population size
duration <- 720 # total number of days
tsteps <- 1 # chunk into 1 day intervals
cs <- .3 # reduction in susceptibility
ci <- .5 # reduction in infectivity
veffs <- (1 - cs) # vaccine efficacy susceptibility
veffi <- (1 - ci) # vaccine efficacy infectivity
beta <- 0.25 # infection rate 1 day^-1
gamma <- 0.1 # recovery rate 1 day^-1
pvacc <- .60 # vaccination coverage aka p
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0) # threshold for vaccine 100% efficacy, version 1
# pvacc <- pvacc_thresh # set vaccination coverage to cvt
peff <- pvacc * veff # effective vaccine coverage
pvacc_thresh <- (1 - (1/R0))/veffs # recalc results to consider efficacy-susceptibility, version 2
pvacc_thresh <- (1 - (1/R0))/(1 - (cs * ci)) # recalc results to consider efficacy-susceptibility & infectivity, version 3
# pvacc <- pvacc_thresh # set vaccination coverage to recalc cvt
(parameters <- c(
peff = peff,
cs = cs,
ci = ci,
veffs = veffs,
veffi = veffi,
beta = beta,
gamma = gamma,
R0 = R0,
pvacc = pvacc,
pvacc_thresh = pvacc_thresh
))## peff cs ci veffs veffi beta
## 0.4200000 0.3000000 0.5000000 0.7000000 0.5000000 0.2500000
## gamma R0 pvacc pvacc_thresh
## 0.1000000 2.5000000 0.6000000 0.7058824initial_state_values <- c(S = round((N-1) * (1-pvacc)),
I = 1,
Iv = 0,
R = 0,
V = round((N-1) * (pvacc)) )
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_c4b <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + Iv + R + V
lambda <- (beta * I/N) + (ci * beta * Iv/N)
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dIv <- (cs * (lambda * V)) -(gamma * Iv)
dR <- (gamma * I) +(gamma * Iv)
dV <- -(cs * (lambda * V))
return(list(c(dS, dI, dIv, dR, dV)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_c4b,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + Iv + R + V),digits=5)*100,
preval_Inf = round(I/(S+ I + Iv + R + V),digits=5)*100,
preval_InfVac = round(Iv/(S+ I + Iv + R + V),digits=5)*100,
propor_Re = round(R/(S + I + Iv + R + V),digits=5)*100,
propor_Vac = round(V/(S + I + Iv + R + V),digits=5)*100,
count_Inf2 = I+ Iv,
preval_Inf2 = round((I+ Iv)/(S+ I + Iv + R + V),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + Iv + R + V)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output %>%
arrange(- (I + Iv), time) %>%
head(1)## time S I Iv R V still_Su preval_Inf
## 1 478 314958.3 8079.541 4224.824 114256.2 558481.1 31.496 0.808
## preval_InfVac propor_Re propor_Vac count_Inf2 preval_Inf2 Reff
## 1 0.422 11.426 55.848 12304.36 1.23 0.7873959print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S I Iv R V still_Su preval_Inf
## 1 1 399999.9 1.001097 0.04331681 0.1022301 599999 40 0
## preval_InfVac propor_Re propor_Vac count_Inf2 preval_Inf2 Reff
## 1 0 0 60 1.044413 0 0.9999997print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I Iv R V still_Su preval_Inf
## 721 720 246230.5 250.938 158.4449 234637.7 518722.5 24.623 0.025
## preval_InfVac propor_Re propor_Vac count_Inf2 preval_Inf2 Reff
## 721 0.016 23.464 51.872 409.3829 0.041 0.6155762print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -preval_InfVac, -propor_Re, -propor_Vac, -Reff, -count_Inf2, -preval_Inf2) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","coral2", "green","steelblue")) +
scale_shape_manual(values = c(0,4,3,1,20)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" cs of", parameters['cs'],
" ci of", parameters['ci'],
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output %>%
select(time, Reff) %>%
filter(time < 120) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" cs of", parameters['cs'],
" ci of", parameters['ci'],
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
Question: For this vaccine, what is the minimum vaccination coverage required to interrupt transmission, i.e. to bring Reff below 1? Calculate then confirm using your model.
”
Relating Reff to R0:
\[\begin{align} R_{eff} & = (1-p) \times R_0 + p c_s c_i \times R_0 \end{align}\]
Solving for pC = p when Reff = 1:
\[\begin{align} p_c & = \frac{1-\frac{1}{R_0}}{1- c_s \times c_i} \\ p_c & = \frac{1-\frac{1}{\frac{0.25}{0.1}}}{1- 0.3 \times 0.5} \\ p_c & = 0.71 \end{align}\]
That is, we need to vaccinate at least 71% of the population in order to interrupt transmission (Reff < 1). ”
- Using critical vaccination threshold calculated with vaccination efficacy for susceptibility and infectivity prevents an epidemic:
nicesubtitle <- "SIR Model v4b2 vaccine w/ combined leaky effects"
print(nicesubtitle)## [1] "SIR Model v4b2 vaccine w/ combined leaky effects"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 1000000 # population size
duration <- 720 # total number of days
tsteps <- 1 # chunk into 1 day intervals
cs <- .3 # reduction in susceptibility
ci <- .5 # reduction in infectivity
veffs <- (1 - cs) # vaccine efficacy susceptibility
veffi <- (1 - ci) # vaccine efficacy infectivity
beta <- 0.25 # infection rate 1 day^-1
gamma <- 0.1 # recovery rate 1 day^-1
pvacc <- .60 # vaccination coverage aka p
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0) # threshold for vaccine 100% efficacy, version 1
# pvacc <- pvacc_thresh # set vaccination coverage to cvt
peff <- pvacc * veff # effective vaccine coverage
pvacc_thresh <- (1 - (1/R0))/veffs # recalc results to consider efficacy-susceptibility, version 2
pvacc_thresh <- (1 - (1/R0))/(1 - (cs * ci)) # recalc results to consider efficacy-susceptibility & infectivity, version 3
pvacc <- pvacc_thresh # set vaccination coverage to recalc cvt
(parameters <- c(
peff = peff,
cs = cs,
ci = ci,
veffs = veffs,
veffi = veffi,
beta = beta,
gamma = gamma,
R0 = R0,
pvacc = pvacc,
pvacc_thresh = pvacc_thresh
))## peff cs ci veffs veffi beta
## 0.4200000 0.3000000 0.5000000 0.7000000 0.5000000 0.2500000
## gamma R0 pvacc pvacc_thresh
## 0.1000000 2.5000000 0.7058824 0.7058824initial_state_values <- c(S = round((N-1) * (1-pvacc)),
I = 1,
Iv = 0,
R = 0,
V = round((N-1) * (pvacc)) )
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_c4b <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + Iv + R + V
lambda <- (beta * I/N) + (ci * beta * Iv/N)
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dIv <- (cs * (lambda * V)) -(gamma * Iv)
dR <- (gamma * I) +(gamma * Iv)
dV <- -(cs * (lambda * V))
return(list(c(dS, dI, dIv, dR, dV)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_c4b,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + Iv + R + V),digits=5)*100,
preval_Inf = round(I/(S+ I + Iv + R + V),digits=5)*100,
preval_InfVac = round(Iv/(S+ I + Iv + R + V),digits=5)*100,
propor_Re = round(R/(S + I + Iv + R + V),digits=5)*100,
propor_Vac = round(V/(S + I + Iv + R + V),digits=5)*100,
count_Inf2 = I+ Iv,
preval_Inf2 = round((I+ Iv)/(S+ I + Iv + R + V),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + Iv + R + V)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output %>%
arrange(- (I + Iv), time) %>%
head(1)## time S I Iv R V still_Su preval_Inf
## 1 92 294110.2 0.7352425 0.5293104 11.3702 705877.1 29.411 0
## preval_InfVac propor_Re propor_Vac count_Inf2 preval_Inf2 Reff
## 1 0 0.001 70.588 1.264553 0 0.7352756print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S I Iv R V still_Su preval_Inf preval_InfVac propor_Re
## 1 0 294117 1 0 0 705882 29.412 0 0 0
## propor_Vac count_Inf2 preval_Inf2 Reff
## 1 70.588 1 0 0.7352925print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I Iv R V still_Su preval_Inf
## 721 720 294064.2 0.7313098 0.5266093 90.62938 705844 29.406 0
## preval_InfVac propor_Re propor_Vac count_Inf2 preval_Inf2 Reff
## 721 0 0.009 70.584 1.257919 0 0.7351604print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -preval_InfVac, -propor_Re, -propor_Vac, -Reff, -count_Inf2, -preval_Inf2) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","coral2", "green","steelblue")) +
scale_shape_manual(values = c(0,4,3,1,20)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" cs of", parameters['cs'],
" ci of", parameters['ci'],
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output %>%
select(time, Reff) %>%
filter(time < 120) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" cs of", parameters['cs'],
" ci of", parameters['ci'],
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
- Using an increased, rounded-up critical vaccination threshold of 0.84 (because last run using calculated cvt had a peak at day 92 with a peak of I + Iv = 1.3, and .84 stops the peak at day 1):
nicesubtitle <- "SIR Model v4b3 vaccine w/ combined leaky effects"
print(nicesubtitle)## [1] "SIR Model v4b3 vaccine w/ combined leaky effects"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 1000000 # population size
duration <- 720 # total number of days
tsteps <- 1 # chunk into 1 day intervals
cs <- .3 # reduction in susceptibility
ci <- .5 # reduction in infectivity
veffs <- (1 - cs) # vaccine efficacy susceptibility
veffi <- (1 - ci) # vaccine efficacy infectivity
beta <- 0.25 # infection rate 1 day^-1
gamma <- 0.1 # recovery rate 1 day^-1
pvacc <- .84 # vaccination coverage aka p
R0 <- beta / gamma
pvacc_thresh <- 1 - (1/R0) # threshold for vaccine 100% efficacy, version 1
# pvacc <- pvacc_thresh # set vaccination coverage to cvt
peff <- pvacc * veff # effective vaccine coverage
pvacc_thresh <- (1 - (1/R0))/veffs # recalc results to consider efficacy-susceptibility, version 2
pvacc_thresh <- (1 - (1/R0))/(1 - (cs * ci)) # recalc results to consider efficacy-susceptibility & infectivity, version 3
# pvacc <- pvacc_thresh # set vaccination coverage to recalc cvt
(parameters <- c(
peff = peff,
cs = cs,
ci = ci,
veffs = veffs,
veffi = veffi,
beta = beta,
gamma = gamma,
R0 = R0,
pvacc = pvacc,
pvacc_thresh = pvacc_thresh
))## peff cs ci veffs veffi beta
## 0.5880000 0.3000000 0.5000000 0.7000000 0.5000000 0.2500000
## gamma R0 pvacc pvacc_thresh
## 0.1000000 2.5000000 0.8400000 0.7058824initial_state_values <- c(S = round((N-1) * (1-pvacc)),
I = 1,
Iv = 0,
R = 0,
V = round((N-1) * (pvacc)) )
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model_c4b <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + Iv + R + V
lambda <- (beta * I/N) + (ci * beta * Iv/N)
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dIv <- (cs * (lambda * V)) -(gamma * Iv)
dR <- (gamma * I) +(gamma * Iv)
dV <- -(cs * (lambda * V))
return(list(c(dS, dI, dIv, dR, dV)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model_c4b,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + Iv + R + V),digits=5)*100,
preval_Inf = round(I/(S+ I + Iv + R + V),digits=5)*100,
preval_InfVac = round(Iv/(S+ I + Iv + R + V),digits=5)*100,
propor_Re = round(R/(S + I + Iv + R + V),digits=5)*100,
propor_Vac = round(V/(S + I + Iv + R + V),digits=5)*100,
count_Inf2 = I+ Iv,
preval_Inf2 = round((I+ Iv)/(S+ I + Iv + R + V),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + Iv + R + V)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output %>%
arrange(- (I + Iv), time) %>%
head(1)## time S I Iv R V still_Su preval_Inf
## 1 1 160000 0.9423567 0.05909123 0.1000981 839998.9 16 0
## preval_InfVac propor_Re propor_Vac count_Inf2 preval_Inf2 Reff
## 1 0 0 84 1.001448 0 0.3999999print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S I Iv R V still_Su preval_Inf preval_InfVac propor_Re
## 1 0 160000 1 0 0 839999 16 0 0 0
## propor_Vac count_Inf2 preval_Inf2 Reff
## 1 84 1 0 0.4print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I Iv R V still_Su
## 721 720 159998.6 7.251942e-10 1.142187e-09 4.613991 839996.8 16
## preval_Inf preval_InfVac propor_Re propor_Vac count_Inf2 preval_Inf2
## 721 0 0 0 84 1.867381e-09 0
## Reff
## 721 0.3999965print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
select(-still_Su, -preval_Inf, -preval_InfVac, -propor_Re, -propor_Vac, -Reff, -count_Inf2, -preval_Inf2) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","coral2", "green","steelblue")) +
scale_shape_manual(values = c(0,4,3,1,20)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" cs of", parameters['cs'],
" ci of", parameters['ci'],
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output %>%
select(time, Reff) %>%
filter(time < 120) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" cs of", parameters['cs'],
" ci of", parameters['ci'],
" p of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
Manual calibration of an SIR model (part 1)
The next step is to make sure that your model is able to capture real-world data. In this exercise you will learn how to do this, with some simple examples. The cell below contains the epidemic data you will be working with, with a column for the time in days, and a column for the prevalence of infection on each day.
Plot the epidemic curve.
Question: Based on the plot, what does the data represent?
- This plot represents the number of infected people (Y axis) over time (x axis).
dataprov <- data.frame(time = 1:5,
number_infected = c(3,8,26,76,225))
dataprov %>%
select(time, number_infected) %>%
ggplot(aes(x=time, y=number_infected)) +
geom_line(linetype="dotdash", color="darkorchid") +
geom_point(color="darkorchid", shape=5, show.legend = FALSE) +
xlab("Time") +
ylab("Number infected, provided") +
labs(
subtitle= "Epidemic Curve, Modelling - Part 1") +
theme(legend.position="bottom")
Calibrate a simple SIR model to this data, i.e. find the combination of \(\beta\) and \(\gamma\) that best reproduces the data. Perform a manual calibration by varying the parameter values manually and compare the model output with the data visually. Through trial-and-error, you find the best parameter values. This manual approach should help to familiarise you with what is happening when you calibrate a model.
Question: Based on the code, what is the total size of the population you are modelling? For plotting, which variables in the model output correspond to the variables in the dataset?
- The provided code snippet provides for N = 763, where we start with one infected individual. However, the provided dataset used as reference starts with 3 infected persons, so this shall be adjusted on the SIR initial conditions, accordingly.
- I (Infected) represents the number infected on the model, which should be tracked against the number_infected variable of the provided datatset.
# INPUT
initial_state_values <- c(S = 762,
I = 1,
R = 0)
times <- seq(from = 0, to = 6, by = 0.1)
# MODEL FUNCTION
sir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S+I+R
lambda <- beta * I/N
# The differential equations
dS <- -lambda * S
dI <- lambda * S - gamma * I
dR <- gamma * I
# Output
return(list(c(dS, dI, dR)))
})Write code to specify values for \(\beta\) and \(\gamma\), solve the model, and plot the model output and the data in the same graph. This plot will show how well your model prediction fits to the data. Try different values for \(\beta\) and \(\gamma\) until you find a combination that produces an output closely matching the observed data.
Question: Which value for \(\beta\) and value for \(\gamma\) did you find to give the best fit to the data?
Assuming an initial prevalence of I0 = 3: using the values of beta \(\beta\) = 1.1 and gamma \(\gamma\) = 0.1 for an R0 R0 of 11, I (Infected) was modeled to be slightly aggressive (higher) than the provided data, but will probably converge on day 6. Otherwise, \(\beta\) = 1.2 and \(\gamma\) = 0.2 will more closely match the provided data through day 5, and undercount by day 6 (guessing a relatively straight line projection).*
Assuming an initial prevalence of I0 = 1: using the values of beta \(\beta\) = 1.16 and gamma \(\gamma\) = 0.01 for an R0 R0 of 116 visually fits relatively well through day 5.
nicesubtitle <- "SIR Model v1a2 SIR Basic Model, manual calibration"
print(nicesubtitle)## [1] "SIR Model v1a2 SIR Basic Model, manual calibration"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 763 # population size
duration <- 6 # total number of days
tsteps <- 0.1 # chunk into 1 day intervals
beta <- 1.1 # infection rate day^-1
gamma <- 0.1 # recovery rate day^-1
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate for untreated I
R0 = R0
)) ## beta gamma R0
## 1.1 0.1 11.0initial_state_values <- c(S = N-3,
I = 3,
R = 0)
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dR <- (gamma * I)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R),digits=5)*100,
preval_Inf = round(I/(S+ I + R),digits=5)*100,
propor_Re = round(R/(S + I + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output %>%
arrange(-I, time) %>%
head(1)## time S I R still_Su preval_Inf propor_Re Reff
## 1 6 289.131 406.8333 67.03569 37.894 53.32 8.786 4.168337print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## [1] time S I R still_Su preval_Inf propor_Re
## [8] Reff
## <0 rows> (or 0-length row.names)print("last record for the run: ")## [1] "last record for the run: "output %>%
arrange(time) %>%
tail(1)## time S I R still_Su preval_Inf propor_Re Reff
## 61 6 289.131 406.8333 67.03569 37.894 53.32 8.786 4.168337print("Plotting the number of people in each compartment over time")## [1] "Plotting the number of people in each compartment over time"output %>%
filter(time <= 5) %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") 
print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output %>%
filter(time <= 5) %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output %>%
filter(time <= 5) %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
# Compare model I to infection numbers provided
output %>%
left_join(dataprov) %>%
select(time, I, number_infected) %>%
melt(id = "time") %>%
filter(!is.na(value)) %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dotdash") +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("red","darkorchid")) +
scale_shape_manual(values = c(4,5)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" I0 of", initial_state_values['I']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") ## Joining, by = "time"
nicesubtitle <- "SIR Model v1a3 SIR Basic Model, manual calibration"
print(nicesubtitle)## [1] "SIR Model v1a3 SIR Basic Model, manual calibration"print("Trial-and-error state values and parameters, patch only")## [1] "Trial-and-error state values and parameters, patch only"# MODEL INPUTS:
beta <- 1.2 # infection rate day^-1
gamma <- 0.2 # recovery rate day^-1
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate for untreated I
R0 = R0
)) ## beta gamma R0
## 1.2 0.2 6.0# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R),digits=5)*100,
preval_Inf = round(I/(S+ I + R),digits=5)*100,
propor_Re = round(R/(S + I + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R)) )
# Compare model I to infection numbers provided
output %>%
left_join(dataprov) %>%
select(time, I, number_infected) %>%
melt(id = "time") %>%
filter(!is.na(value)) %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dotdash") +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("red","darkorchid")) +
scale_shape_manual(values = c(4,5)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" I0 of", initial_state_values['I']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") ## Joining, by = "time"
nicesubtitle <- "SIR Model v1a4 SIR Basic Model, manual calibration"
print(nicesubtitle)## [1] "SIR Model v1a4 SIR Basic Model, manual calibration"print("Trial-and-error state values and parameters, patch only")## [1] "Trial-and-error state values and parameters, patch only"# MODEL INPUTS:
beta <- 0.96 # infection rate day^-1
gamma <- 0.03 # recovery rate day^-1
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate for untreated I
R0 = R0
)) ## beta gamma R0
## 0.96 0.03 32.00# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R),digits=5)*100,
preval_Inf = round(I/(S+ I + R),digits=5)*100,
propor_Re = round(R/(S + I + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R)) )
# Compare model I to infection numbers provided
output %>%
left_join(dataprov) %>%
select(time, I, number_infected) %>%
melt(id = "time") %>%
filter(!is.na(value)) %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dotdash") +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("red","darkorchid")) +
scale_shape_manual(values = c(4,5)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" I0 of", initial_state_values['I']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") ## Joining, by = "time"
- Settle back to initial I I0 = 1 to get a different set of values for next exercise.
- Using the values of beta \(\beta\) = 1.16 and gamma \(\gamma\) = 0.01.
nicesubtitle <- "SIR Model v1a5 SIR Basic Model, manual calibration, diff Init Values"
print(nicesubtitle)## [1] "SIR Model v1a5 SIR Basic Model, manual calibration, diff Init Values"print("Trial-and-error state values and parameters, patch only")## [1] "Trial-and-error state values and parameters, patch only"initial_state_values <- c(S = N-1,
I = 1,
R = 0)
# MODEL INPUTS:
beta <- 1.16 # infection rate day^-1
gamma <- 0.01 # recovery rate day^-1
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate for untreated I
R0 = R0
)) ## beta gamma R0
## 1.16 0.01 116.00# MODEL OUTPUT (solving the differential equations):
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R),digits=5)*100,
preval_Inf = round(I/(S+ I + R),digits=5)*100,
propor_Re = round(R/(S + I + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R)) )
# Compare model I to infection numbers provided
output %>%
left_join(dataprov) %>%
select(time, I, number_infected) %>%
melt(id = "time") %>%
filter(!is.na(value)) %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dotdash") +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("red","darkorchid")) +
scale_shape_manual(values = c(4,5)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" I0 of", initial_state_values['I']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") ## Joining, by = "time"
Manual calibration of an SIR model (part 2)
Most people found different combinations of \(\beta\) and \(\gamma\) in the last activity, because there is not enough data in the early epidemic curve to identify unique values for the two parameters separately. The particular outbreak we are looking at is actually already over, and we have data for the following days.
Here is the data for the full epidemic (not just the initial growth), saved in the “data” variable, with the number of infected people recorded over 14 days:
dataprov2 <- data.frame(time = 1:14,
number_infected = c(3,8,26,76,225,298,258,233,189,128,68,29,14,4))
dataprov2 %>%
select(time, number_infected) %>%
ggplot(aes(x=time, y=number_infected)) +
geom_line(linetype="dotdash", color="darkorchid") +
geom_point(color="darkorchid", shape=5, show.legend = FALSE) +
xlab("Time") +
ylab("Number infected, provided") +
labs(
subtitle= "Epidemic Curve, Modelling - Part 2") +
theme(legend.position="bottom")
Now, let’s see what happens if we run the model again, with the parameter values we found before, but for a longer time (you can fill in your own values below, if you like). Does it still also fit the later part of the epidemic curve?
Reuse parameters from most recent run, but use longer duration:
Initial I I0 = 1
Beta \(\beta\) = 1.16
gamma \(\gamma\) = 0.01
The prior parameters no longer fit the rest of the (new, after day 5) data well.
nicesubtitle <- "SIR Model v1a6 SIR Basic Model, manual calibration"
print(nicesubtitle)## [1] "SIR Model v1a6 SIR Basic Model, manual calibration"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 763 # population size
duration <- 14 # total number of days
tsteps <- 0.2 # chunk into 1 day intervals
beta <- 1.16 # infection rate day^-1
gamma <- 0.01 # recovery rate day^-1
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate for untreated I
R0 = R0
)) ## beta gamma R0
## 1.16 0.01 116.00initial_state_values <- c(S = N-1,
I = 1,
R = 0)
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dR <- (gamma * I)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
output2 <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R),digits=5)*100,
preval_Inf = round(I/(S+ I + R),digits=5)*100,
propor_Re = round(R/(S + I + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output2 %>%
arrange(-I, time) %>%
head(1)## time S I R still_Su preval_Inf propor_Re Reff
## 1 10 6.286704 725.1572 31.55605 0.824 95.04 4.136 0.9557768print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output2 %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S I R still_Su preval_Inf propor_Re Reff
## 1 10 6.286704 725.1572 31.55605 0.824 95.04 4.136 0.9557768print("last record for the run: ")## [1] "last record for the run: "output2 %>%
arrange(time) %>%
tail(1)## time S I R still_Su preval_Inf propor_Re Reff
## 71 14 0.08098028 702.7374 60.18159 0.011 92.102 7.887 0.01231155print("Plotting the number of people in each compartment over time")## [1] "Plotting the number of people in each compartment over time"output2 %>%
filter(time <= 5) %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") 
print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output2 %>%
filter(time <= 5) %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output2 %>%
filter(time <= 5) %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
# Compare model I to infection numbers provided
output2 %>%
left_join(dataprov2) %>%
select(time, I, number_infected) %>%
melt(id = "time") %>%
filter(!is.na(value)) %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dotdash") +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("red","darkorchid")) +
scale_shape_manual(values = c(4,5)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" I0 of", initial_state_values['I']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") ## Joining, by = "time"
The model completely overestimates the number of infected people at later timepoints. The last datapoint the model was calibrated to was already close to the peak of the epidemic, but was nothing in the earlier, limited data to suggest this. Limited, early data can give rise to a range of projected epidemic curves, not all of them consistent with how the epidemic actually developed.
Find values of \(\beta\) and \(\gamma\) that match the model as closely as possible to this new data, for the full epidemic curve? Calibrate the model to the full dataset provided here.
Question: Which value for \(\beta\) and value for \(\gamma\) did you find to give the best fit to the data?
Visually, the best fit by trial-and-error was:
- Initial I I0 = 1
- Beta \(\beta\) = 1.72
- gamma \(\gamma\) = 0.43
Background information on the dataset:
The dataset we use here is from a real influenza outbreak, which occurred in an English boarding school in 1978. The number infected actually represents the number of boys that were confined to bed each day.
The numbers presented here were taken from De Vries et al. (2006) and if you are interested, you can find out more about the outbreak from the original publication in the British Medical Journal.
References:
Anonymous. 1978. Influenza in a boarding school. British Medical Journal 1:578.
G. De Vries, T. Hillen, M. Lewis, J. Mueller, and B. Schoenfisch. 2006. A Course in Mathematical Biology: Quantitative Modeling with Mathematical and Computational Methods. Society for Industrial and Applied Mathematics.
nicesubtitle <- "SIR Model v1a7 SIR Basic Model, manual calibration"
print(nicesubtitle)## [1] "SIR Model v1a7 SIR Basic Model, manual calibration"print("Trial-and-error state values and parameters, patch only")## [1] "Trial-and-error state values and parameters, patch only"# MODEL INPUTS:
beta <- 1.72 # infection rate day^-1
gamma <- 0.43 # recovery rate day^-1
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate for untreated I
R0 = R0
)) ## beta gamma R0
## 1.72 0.43 4.00# MODEL OUTPUT (solving the differential equations):
output2 <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R),digits=5)*100,
preval_Inf = round(I/(S+ I + R),digits=5)*100,
propor_Re = round(R/(S + I + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R)) )
# Compare model I to infection numbers provided
output2 %>%
left_join(dataprov2) %>%
select(time, I, number_infected) %>%
melt(id = "time") %>%
filter(!is.na(value)) %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dotdash") +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("red","darkorchid")) +
scale_shape_manual(values = c(4,5)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" I0 of", initial_state_values['I']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") ## Joining, by = "time"
- Provided solution uses \(\beta\) = 1.70 and \(\gamma\) = 0.45
nicesubtitle <- "SIR Model v1a8 SIR Basic Model, manual calibration"
print(nicesubtitle)## [1] "SIR Model v1a8 SIR Basic Model, manual calibration"print("Trial-and-error state values and parameters, patch only")## [1] "Trial-and-error state values and parameters, patch only"# MODEL INPUTS:
beta <- 1.70 # infection rate day^-1
gamma <- 0.45 # recovery rate day^-1
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate for untreated I
R0 = R0
)) ## beta gamma R0
## 1.700000 0.450000 3.777778# MODEL OUTPUT (solving the differential equations):
output2 <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R),digits=5)*100,
preval_Inf = round(I/(S+ I + R),digits=5)*100,
propor_Re = round(R/(S + I + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R)) )
# Compare model I to infection numbers provided
output2 %>%
left_join(dataprov2) %>%
select(time, I, number_infected) %>%
melt(id = "time") %>%
filter(!is.na(value)) %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dotdash") +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("red","darkorchid")) +
scale_shape_manual(values = c(4,5)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3),
" I0 of", initial_state_values['I']),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") ## Joining, by = "time"
Writing a sum-of-squares function in R, for SIR model fits
Introduction
When you are trying to fit a model to data, you use a distance function to assess, quantitatively how “far” the model is from the data. By calculating the output of this function for a range of parameter sets, you can find which of those parameter sets generates the model which is the closest fit to the data.
You can calculate the value of the distance function by hand for each model you try - but because this means doing the same calculation repeatedly with different inputs, this is a good task for which using a function would be efficient.
So here, you will create a function yourself, to
calculate and return the sum-of-squares value (SSQ) for the fit of a simple SIR model, parameterised with given values of \(\beta\) and \(\gamma\), to any dataset.
use this function to find the sum-of-squares of models fit to a simple datset (the flu dataset you may have used in other activities), with the following parameters:
- \(\beta\) = 1.15, \(\gamma\) = 0.02
- \(\beta\) = 1.7, \(\gamma\) = 0.45
- \(\beta\) = 1.9, \(\gamma\) = 0.6
- \(\beta\) = 1.15, \(\gamma\) = 0.02
Put very simply, your function will have this structure:
SIR_SSQ <- function(arguments) { ### fill in your arguments
# Calculate model output
### YOUR CODE HERE ###
# Calculate sum-of-squares (SSQ) of model fit
### YOUR CODE HERE ###
return(SSQ)
}Calculating the SSQ gives a quantitative measure of the distance from your model output to the data.
Last week’s activities demonstrate manual calibration of a model, visually checking the model ouputs from a small number of parameter combinations. It would quickly become infeasible to test enough combinations of parameters in this way, to find the combination producing the best available model fit. So automating the process of simulating each model and calculating the SSQ, giving a quantitative distance measure instead of a visual check for each model, makes this more efficient.
So in this activity we will create a function which can: simulate the model for a given combination of parameters, compare the output to data and calculate the SSQ.
Later in this course you will learn how to use optimisation functions to take this a step further and automatically find the combination of parameters giving the best fit, as defined by your chosen distance function. One of these functions, optim(), is introduced in the activity following this one in this module.
With this in mind, think about
* what else you will need to put within the body of the function, apart from the sum-of-squares calculation
* what inputs you will give the function and how to arrange them into arguments.
Tips:
When using
ode(), you have an example of a function as an argument of another function. The function that you feed intoode()gets the values for its arguments from theode()function.Keep track of the names of your arguments, so you can see that they match with what the “inner” function is expecting.
To make your function applicable to any dataset and combination of parameters you choose, the function will need to read in the data as one of its arguments.
To calculate the SSQ, you will have to calculate the difference between model output and data, at only the timepoints for which data is available. How can you achieve this?
Task:
create a function which calculates and returns the sum-of-squares value (SSQ) for the fit of a simple SIR model, parameterised with given values of \(\beta\) and \(\gamma\), to any dataset.
use this function to find the sum-of-squares of models fit to a simple datset (the flu dataset you may have used in other activities), with the following parameters:
- \(\beta\) = 1.15, \(\gamma\) = 0.02 ~ ssq 2507764.33
- \(\beta\) = 1.7, \(\gamma\) = 0.45 ~ ssq 4630.30
- \(\beta\) = 1.9, \(\gamma\) = 0.6 ~ ssq 34243.68
- \(\beta\) = 1.15, \(\gamma\) = 0.02 ~ ssq 2507764.33
# SIR MODEL FUNCTION
sir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dR <- (gamma * I)
return(list(c(dS, dI, dR)))
})
}
# SSQ FUNCTION (assumes sir_model model already loaded)
# parameters - vector with named elements for sir_model
# idat - df or list containing vectors number_infected (I) and time
SIR_SSQ <- function(parameters, idat) {
# MODEL OUTPUT (solving the differential equations):
oderesult <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters))
# sum-of-squares (SSQ) of model fit requires idat
idat <- na.omit(idat)
deltas2 <- (oderesult$I[oderesult$time %in% idat$time] - idat$number_infected)^2
SSQ <- sum(deltas2)
return(SSQ)
}
## load data
flu_data <- data.frame(time = 1:14,
number_infected = c(3,8,26,76,225,298,258,233,189,128,68,29,14,4))
initial_state_values <- c(S = 762, I = 1, R = 0)
## calculate the sum-of-squares for an SIR model fit to these data where:
### beta = 1.15, gamma = 0.02
### beta = 1.7, gamma = 0.45
### beta = 1.9, gamma = 0.6
# populate time for 14 days, dense
times <- seq(from = 0, to = 14, by = 0.1)
# test 1 with beta and gamma in a vector 'parms'
parms <- c(beta = 1.15, gamma = 0.02)
ssq1 <- SIR_SSQ(parameters = parms,
idat = flu_data)
print(paste0("ssq = ", ssq1, ", where beta = ", parms['beta'], ", gamma = ", parms['gamma']))## [1] "ssq = 2507764.33057048, where beta = 1.15, gamma = 0.02"# test 2 with a beta and gamma entered in function
ssq2 <- SIR_SSQ(parameters =c(beta =1.7, gamma =0.45),
idat =flu_data)
print(paste0("ssq = ", ssq2, ", where beta = 1.7, gamma = 0.45"))## [1] "ssq = 4630.30225929158, where beta = 1.7, gamma = 0.45"# test 3 with beta and gamma in a vector 'parms'
parms <- c(beta = 1.9, gamma = 0.6)
ssq3 <- SIR_SSQ(parameters = parms,
idat = flu_data)
print(paste0("ssq = ", ssq3, ", where beta = ", parms['beta'], ", gamma = ", parms['gamma']))## [1] "ssq = 34243.6805666932, where beta = 1.9, gamma = 0.6"Further work:
Modify your SSQ function so it would take in any type of model equation into the ode() function:
# SIR MODEL FUNCTION
sir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dR <- (gamma * I)
return(list(c(dS, dI, dR)))
})
}
# SSQ FUNCTION (assumes sir_model model already loaded)
# parameters - vector with named elements for sir_model
# idat - df or list containing vectors number_infected (I) and time
SIR_SSQv2 <- function(func, parameters, idat) {
# MODEL OUTPUT (solving the differential equations):
oderesult <- as.data.frame(ode(y = initial_state_values,
times = times,
func = func,
parms = parameters))
# sum-of-squares (SSQ) of model fit requires idat
idat <- na.omit(idat)
deltas2 <- (oderesult$I[oderesult$time %in% idat$time] - idat$number_infected)^2
SSQ <- sum(deltas2)
return(SSQ)
}
## load data
flu_data <- data.frame(time = 1:14,
number_infected = c(3,8,26,76,225,298,258,233,189,128,68,29,14,4))
initial_state_values <- c(S = 762, I = 1, R = 0)
## calculate the sum-of-squares for an SIR model fit to these data where:
### beta = 1.15, gamma = 0.02
### beta = 1.7, gamma = 0.45
### beta = 1.9, gamma = 0.6
# populate time for 14 days, dense
times <- seq(from = 0, to = 14, by = 0.1)
# test 1 with beta and gamma in a vector 'parms'
parms <- c(beta = 1.15, gamma = 0.02)
ssq1 <- SIR_SSQv2(func = sir_model,
parameters = parms,
idat = flu_data)
print(paste0("ssq = ", ssq1, ", where beta = ", parms['beta'], ", gamma = ", parms['gamma']))## [1] "ssq = 2507764.33057048, where beta = 1.15, gamma = 0.02"# test 2 with a beta and gamma entered in function
ssq2 <- SIR_SSQv2(func = sir_model,
parameters =c(beta =1.7, gamma =0.45),
idat =flu_data)
print(paste0("ssq = ", ssq2, ", where beta = 1.7, gamma = 0.45"))## [1] "ssq = 4630.30225929158, where beta = 1.7, gamma = 0.45"# test 3 with beta and gamma in a vector 'parms'
parms <- c(beta = 1.9, gamma = 0.6)
ssq3 <- SIR_SSQv2(func = sir_model,
parameters = parms,
idat = flu_data)
print(paste0("ssq = ", ssq3, ", where beta = ", parms['beta'], ", gamma = ", parms['gamma']))## [1] "ssq = 34243.6805666932, where beta = 1.9, gamma = 0.6"Automated least-squares calibration: Using optim(), R’s built-in optimisation function
Introduction
You should now have a working function which takes in data and parameters for an SIR model and calculates the sum-of-squares for the model fit. This is now a function which you can put into optim() to find the \(\beta\) and \(\gamma\) for the best SIR model fit to a dataset, as measured by least squares.
Using optim():
Arguments it needs at a minimum are:
par : initial values for the parameters to be optimised over
fn : a function to be minimised
So optim() takes a function (which you provide in the fn argument) for which it needs to find the minimum (or maximum) value. The output value of the function depends on the values of the parameters, which optim() will change as it works through the optimisation process; you need to give it some starting values in the par argument.
What other argument might you need to give?
- Starting values for parameters beta and gamma
Tasks
- Use
optim()to find the beta and gamma for the best SIR model to fit an example dataset
- Plot the example dataset against your “best fit” model
# SIR MODEL FUNCTION
sir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dR <- (gamma * I)
return(list(c(dS, dI, dR)))
})
}
# SSQ FUNCTION (assumes sir_model model already loaded)
# parameters - vector with named elements for sir_model
# idat - df or list containing vectors number_infected (I) and time
SIR_SSQv2 <- function(func, parameters, idat) {
# MODEL OUTPUT (solving the differential equations):
oderesult <- as.data.frame(ode(y = initial_state_values,
times = times,
func = func,
parms = parameters))
# sum-of-squares (SSQ) of model fit requires idat
idat <- na.omit(idat)
deltas2 <- (oderesult$I[oderesult$time %in% idat$time] - idat$number_infected)^2
SSQ <- sum(deltas2)
return(SSQ)
}
## load data
flu_data <- data.frame(time = 1:14,
number_infected = c(3,8,26,76,225,298,258,233,189,128,68,29,14,4))
initial_state_values <- c(S = 762, I = 1, R = 0)
## calculate the sum-of-squares for an SIR model fit to these data where:
### beta = 1.15, gamma = 0.02
# populate time for 14 days, dense
times <- seq(from = 0, to = 14, by = 0.1)
# test 1 with beta and gamma in a vector 'parms'
parms <- c(beta = 1.15, gamma = 0.02)
ssq1 <- SIR_SSQv2(func = sir_model,
parameters = parms,
idat = flu_data)
print(paste0("ssq = ", ssq1, ", where beta = ", parms['beta'], ", gamma = ", parms['gamma']))## [1] "ssq = 2507764.33057048, where beta = 1.15, gamma = 0.02"# choose values to start your optimisation
beta_start <- 1
gamma_start <- 0.5
# run optim
(optimised <- optim(par = c(beta = beta_start
, gamma = gamma_start)
, fn = SIR_SSQv2
, func = sir_model
, idat = flu_data
))## $par
## beta gamma
## 1.6692280 0.4434254
##
## $value
## [1] 4121.945
##
## $counts
## function gradient
## 67 NA
##
## $convergence
## [1] 0
##
## $message
## NULL# examine optim() output and plot "best" model against example dataset
opt_mod <- as.data.frame(ode(y = initial_state_values # named vector of initial state values
, times = times # vector of times
, func = sir_model # your predefined SIR function
, parms = optimised$par
))
nicesubtitle <- "Automated Least-Squares Calibration: Using optim() c2w3_2"
print(nicesubtitle)## [1] "Automated Least-Squares Calibration: Using optim() c2w3_2"# Compare model I to infection numbers provided
opt_mod %>%
left_join(flu_data) %>%
select(time, I, number_infected) %>%
melt(id = "time") %>%
filter(!is.na(value)) %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dotdash") +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("red","darkorchid")) +
scale_shape_manual(values = c(4,5)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", round(optimised$par[1],3),
" gamma of", round(optimised$par[2],3),
" SSQ", round(optimised$value,0),
" R0 of", round(optimised$par[1]/optimised$par[2],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") ## Joining, by = "time"
How calibrations inform policy
In many practical applications of modelling, we don’t know the values of all the parameters. To determine the impact of interventions on an epidemic using mathematical modelling and to inform policy, we first need to estimate the parameter values from the available data, such as the prevalence of infection. Then, we can use the inferred parameter values to make projections about the effect of different interventions.
This exercise will show you how the least-squares approach to calibrating a model can be used in the model development process with the aim of informing policy-relevant questions.
Running the cell below will load a new dataset of the prevalence of infection for an outbreak that lasted 200 days, in a population of 500 people. Calibrate a SIR model to this dataset using the automated least-squares algorithm, and then to use this information to draw conclusions about the amount of vaccination that will be needed, to protect other populations from this disease.
What are the best-fitting values for \(\beta\) and \(\gamma\)?
- beta of 0.16 gamma of 0.05
dataprov3 <- read.csv("m3_nb3_data.csv")
dataprov3 %>%
select(time, number_infected) %>%
ggplot(aes(x=time, y=number_infected)) +
geom_line(linetype="dotdash", color="darkorchid") +
geom_point(color="darkorchid", shape=5, show.legend = FALSE) +
xlab("Time") +
ylab("Number infected, provided") +
labs(
subtitle= "Epidemic Curve, Modelling - c2w3_3") +
theme(legend.position="bottom")
# SIR MODEL FUNCTION
sir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dR <- (gamma * I)
return(list(c(dS, dI, dR)))
})
}
# SSQ FUNCTION (assumes sir_model model already loaded)
# parameters - vector with named elements for sir_model
# idat - df or list containing vectors number_infected (I) and time
SIR_SSQv2 <- function(func, parameters, idat) {
# MODEL OUTPUT (solving the differential equations):
oderesult <- as.data.frame(ode(y = initial_state_values,
times = times,
func = func,
parms = parameters))
# sum-of-squares (SSQ) of model fit requires idat
idat <- na.omit(idat)
deltas2 <- (oderesult$I[oderesult$time %in% idat$time] - idat$number_infected)^2
SSQ <- sum(deltas2)
return(SSQ)
}
# initial_state_values and times
initial_state_values <- c(S = 499, I = 1, R = 0)
times <- seq(from = 0, to = 200, by = 1)
# choose values to start your optimisation
beta_start <- 0.9
gamma_start <- 0.1
# run optim
(optimised <- optim(par = c(beta = beta_start
, gamma = gamma_start)
, fn = SIR_SSQv2
, func = sir_model
, idat = dataprov3
))## $par
## beta gamma
## 0.16035940 0.04973196
##
## $value
## [1] 5160.721
##
## $counts
## function gradient
## 89 NA
##
## $convergence
## [1] 0
##
## $message
## NULL# examine optim() output and plot "best" model against example dataset
opt_mod <- as.data.frame(ode(y = initial_state_values # named vector of initial state values
, times = times # vector of times
, func = sir_model # your predefined SIR function
, parms = optimised$par
))
nicesubtitle <- "Automated Least-Squares Calibration: Using optim() c2w3_3"
print(nicesubtitle)## [1] "Automated Least-Squares Calibration: Using optim() c2w3_3"# Compare model I to infection numbers provided
opt_mod %>%
left_join(dataprov3) %>%
select(time, I, number_infected) %>%
melt(id = "time") %>%
filter(!is.na(value)) %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dotdash") +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("red","darkorchid")) +
scale_shape_manual(values = c(4,5)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", round(optimised$par[1],3),
" gamma of", round(optimised$par[2],3),
" SSQ", round(optimised$value,0),
" R0 of", round(optimised$par[1]/optimised$par[2],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") ## Joining, by = "time"
Based on your estimates parameter values, what would be the critical vaccination threshold required to prevent this epidemic, assuming an all-or-nothing vaccine with 75% efficacy?
” …calculate the critical effective vaccination coverage as:
\[\begin{align} p_{eff} & = 1-\frac{1}{R_0} \\ & = 1-\frac{1}{\frac{\beta}{\gamma}} \\ & = 1-\frac{0.05}{0.16} \\ & = 0.69 \end{align}\]
Since the all-or-nothing vaccine is only 75% effective, we can calculate the critical vaccination threshold as: \[\begin{align} p_{c} & = \frac{p_{eff}}{v_{eff}} \\ & = \frac{0.69}{0.75} \\ & = 0.92 \end{align}\]
To interrupt transmission and bring R0 below 1, 92% of the population would have to be vaccinated. ”
nicesubtitle <- "SIR Model v1a9 SIR Basic Model, check pvacc = pc of .92 \n all-or-nothing vaccine with 75% efficacy"
print(nicesubtitle)## [1] "SIR Model v1a9 SIR Basic Model, check pvacc = pc of .92 \n all-or-nothing vaccine with 75% efficacy"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 500 # population size
duration <- 200 # total number of days
tsteps <- 1 # chunk in days
beta <- 0.16 # infection rate day^-1
gamma <- 0.05 # recovery rate day^-1
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate for untreated I
R0 = R0
)) ## beta gamma R0
## 0.16 0.05 3.20# initial_state_values and times
initial_state_values <- c(S = (N-1) * .08, # 8% are susceptible
I = 1,
R = (N-1) * .92) # 92% are vaccinated / immune
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dR <- (gamma * I)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
output2 <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R),digits=5)*100,
preval_Inf = round(I/(S+ I + R),digits=5)*100,
propor_Re = round(R/(S + I + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R)) )
print("peak infection day when I is at its max: ")## [1] "peak infection day when I is at its max: "output2 %>%
arrange(-I, time) %>%
head(1)## time S I R still_Su preval_Inf propor_Re Reff
## 1 0 39.92 1 459.08 7.984 0.2 91.816 0.255488print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output2 %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S I R still_Su preval_Inf propor_Re Reff
## 1 0 39.92 1 459.08 7.984 0.2 91.816 0.255488print("last record for the run: ")## [1] "last record for the run: "output2 %>%
arrange(time) %>%
tail(1)## time S I R still_Su preval_Inf propor_Re Reff
## 201 200 39.579 0.0005733361 460.4204 7.916 0 92.084 0.2533056print("Plotting the number of people in each compartment over time")## [1] "Plotting the number of people in each compartment over time"output2 %>%
filter(time <= 5) %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") 
print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output2 %>%
filter(time <= 5) %>%
select(-still_Su, -preval_Inf, -propor_Re, -Reff) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line() +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","red","green")) +
scale_shape_manual(values = c(0,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
output2 %>%
filter(time <= 5) %>%
select(time, Reff) %>%
ggplot(aes(x=time, y=Reff)) +
geom_line(linetype="dotdash", color="orange") +
geom_jitter(color="orange", shape=6, show.legend = FALSE) +
xlab("Time (days)") +
ylab("R Effective") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
Calibrating the SIR model using likelihood as a measure of the divergence
We are building up to calibrating the SIR model to our flu outbreak data using likelihood as a measure of the divergence between the model projections and the data.
Even though we are looking at the same outbreak, the dataset only shows the reported cases, and we know that 60% of flu cases are reported.
Load the model function, inputs and datasets. Plot of the dataset, then simulate the SIR model using the parameter values found in the prior manual calibration to the full dataset of the total number infected.
# 1 of 4
# Load the flu dataset of reported cases
dataprov4 <- read.csv("idm2_sir_reported_data.csv")
dataprov4 %>%
select(time, number_reported) %>%
ggplot(aes(x=time, y=number_reported)) +
geom_line(linetype="dotdash", color="cadetblue4") +
geom_point(color="chartreuse3", shape=10, show.legend = FALSE) +
xlab("Time (days)") +
ylab("Number of Reported Cases") +
labs(
subtitle= "Calibrating SIR model using likelihood c2w4_1") +
theme(legend.position="bottom")
nicesubtitle <- "SIR Model v1a10 SIR Basic Model, reported data with provided beta and gamma c2w4_1"
print(nicesubtitle)## [1] "SIR Model v1a10 SIR Basic Model, reported data with provided beta and gamma c2w4_1"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 763 # population size
duration <- 14 # total number of days
tsteps <- 0.1 # chunk in days
beta <- 1.7 # infection rate day^-1
gamma <- 0.45 # recovery rate day^-1
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate
R0 = R0
)) ## beta gamma R0
## 1.700000 0.450000 3.777778# initial_state_values and times
initial_state_values <- c(S = (N-1),
I = 1,
R = 0)
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SIR MODEL FUNCTION
sir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + I + R
lambda <- beta * I/N
dS <- -(lambda * S)
dI <- (lambda * S) -(gamma * I)
dR <- (gamma * I)
return(list(c(dS, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
output2 <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R),digits=5)*100,
preval_Inf = round(I/(S+ I + R),digits=5)*100,
propor_Re = round(R/(S + I + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R)) )
output2 %>%
left_join(dataprov4) %>%
select(time, I, number_reported) %>%
melt(id = "time") %>%
filter(!is.na(value)) %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dotdash") +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("red","cadetblue4")) +
scale_shape_manual(values = c(4,10)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", parameters['gamma'],
" R0 of", round(parameters['R0'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom") ## Joining, by = "time"
We want to calculate the likelihood of the model with these specific parameter values, i.e. the probability of observing these numbers of reported cases given our simulated numbers of infected people.
Question: Calculate the Poisson log-likelihood for the epidemic curve?
# Code block 4:
(LL <- sum(dpois(x = dataprov4$number_reported, lambda = 0.6 * output2$I[output2$time %in% dataprov4$time], log = TRUE)))## [1] -61.51812plot(dpois(x = dataprov4$number_reported, lambda = 0.6 * output2$I[output2$time %in% dataprov4$time], log = TRUE))
The dpois() command calculates the Poisson density - the Poisson distribution is a common choice to model the number of events (e.g. the number of reported infections) in a defined time interval.
For each timepoint, the dpois() command calculates the probability of observing the number of reported cases in the dataset (the x argument) given the number predicted by the model with the specified set of parameters (the lambda argument, not to be confused with the force of infection) – it calculates the likelihood of these parameter values for the given datapoint.
The dataset contains the number of reported cases at each timepoint, but the model simulates the total number of infections. To take this into account when calculating how close the simulation is to the observed data, we need to calculate the number of reported cases the model would predict with our knowledge of the reporting rate, i.e. we need to multiply the simulated total number of infections by 0.6.
We usually want to calculate the log-likelihood, not the likelihood itself. Accordingly, setting “log = TRUE” in dpois() calculates the log-Poisson density.
Finally, as we are working on the log-scale, calculating the overall likelihood for the whole epidemic curve requires summing the individual log-likelihoods using the sum() command.
Performing maximum likelihood estimation
Use the code for calculating the Poisson log likelihood for a given parameter combination, together with the calibration function, to fit the model to our flu outbreak data – by finding the parameter values \(\beta\) and \(\gamma\) that maximise the likelihood. This will be calibrating the model using a likelihood approach.
Change your code to calculate the log-likelihood distance function instead of the sum-of-squares. The optimisation process is similar, except that we need to maximise the likelihood (versus minimise the distance as represented by sum of squares).
Define a function (the distance function described to in the lecture) that simulates the model for a given combination of parameters and calculates the Poisson log-likelihood for the epidemic curve of reported cases.
# 2 of 4
# DISTANCE FUNCTION
loglik_function <- function(parameters, idat) { # param values and dataset
beta <- parameters[1] # first value in params is beta
gamma <- parameters[2] # second value in params is gamma
# Simulate the model with initial conditions and timesteps
oderesult2 <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = c(beta = beta, # beta fr params of loglik_func
gamma = gamma))) # gamma from parameters of loglik_func
# Calculate log-likelihood, accounting for the reporting rate of 60%:
LL <- sum(dpois(x = idat$number_reported, lambda = 0.6 * oderesult2$I[oderesult2$time %in% idat$time], log = TRUE))
return(LL)
}Optimise the function using the optim() command to find the values for beta and gamma giving the highest log-likelihood value as output:
# 3 of 4
# OPTIMISATION:
# optim(par = c(parameters['beta'],parameters['gamma']), # doesn't quite get exact same answer
# optim(par = c(1.7, 0.45),
# optim(par = c(1.0, 0.01),
(optimised2 <- optim(par = c(1.7, 0.1), # starting values for beta and gamma - you should get the same result no matter what but not true
fn = loglik_function,
idat = dataprov4,
control = list(fnscale=-1)) # look for max
)## $par
## [1] 1.6915096 0.4764014
##
## $value
## [1] -59.23995
##
## $counts
## function gradient
## 51 NA
##
## $convergence
## [1] 0
##
## $message
## NULLSince we have the full dataset of the number of infected people, confirm that these parameter values indeed produce a good visual fit to the real data by plotting the model simulation alongside the data, in the cell below.
# 4 of 4
nicesubtitle <- "Automated Maximum Likelihood Calibration: Using optim() c2w4_2"
print(nicesubtitle)## [1] "Automated Maximum Likelihood Calibration: Using optim() c2w4_2"# Load the flu dataset of the total number infected
dataprov5 <- read.csv("idm2_sir_data.csv")
# plot total number infected versus prior reported infected
dataprov5 %>%
left_join(dataprov4) %>%
select(time, number_infected, number_reported) %>%
melt(id = "time") %>%
filter(!is.na(value)) %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dotdash",show.legend = FALSE) +
geom_point() +
scale_color_manual(values = c("cadetblue4","darkorchid")) +
scale_shape_manual(values = c(10,5)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(subtitle= "Reported vs. Full Dataset of Infected - c2w4_2") +
theme(legend.position="bottom") +
theme(legend.title = element_blank()) ## Joining, by = "time"
# simulate model with the best-fitting (max-likelihood) parameters from prior steps
(parameters <- c(beta = optimised2$par[1],
gamma = optimised2$par[2],
R0 = optimised2$par[1]/optimised2$par[2]
))## beta gamma R0
## 1.6915096 0.4764014 3.5505976# MODEL OUTPUT (solving the differential equations):
output2 <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + I + R),digits=5)*100,
preval_Inf = round(I/(S+ I + R),digits=5)*100,
propor_Re = round(R/(S + I + R),digits=5)*100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + I + R)) )
output2 %>%
left_join(dataprov5) %>%
left_join(dataprov4) %>%
select(time, I, number_infected, number_reported) %>%
melt(id = "time") %>%
filter(!is.na(value)) %>%
ggplot(aes(x=time, y=value, color=variable, shape=variable)) +
geom_line(linetype="dotdash",show.legend = FALSE) +
geom_point() +
scale_color_manual(values = c("red", "cadetblue4", "darkorchid")) +
scale_shape_manual(values = c(4, 10, 5)) +
xlab("Time (days)") +
ylab("Number of people") +
labs(title= paste0("Model *I*nfected vs. Reported vs. Full Dataset of Infected",
"\n beta of ", round(parameters['beta'],3),
" gamma of ", round(parameters['gamma'],3),
" R0 of ", round(parameters['R0'],3)),
subtitle= nicesubtitle) +
theme(legend.position="bottom") +
theme(legend.title = element_blank()) ## Joining, by = "time"
## Joining, by = "time"
Assessment
Situation:
A new respiratory virus has begun circulating in another country. So far, epidemiological investigations have yielded the data below. What is an estimate for the average infection rate (the ‘beta’ parameter) for this disease? If a vaccine with 80% effectiveness becomes available before an epidemic begins in your own country, what is the minimum coverage of this vaccine, in order to prevent an outbreak from occurring?
Epidemiological data
- All infected individuals go through an incubation period, which lasts on average for 4 days. During this time, individuals are not infectious, nor do they have any excess mortality risk.
- All infected individuals eventually develop symptoms, and the mean duration of symptoms before recovery or death is 5 days. Symptomatic individuals are infectious, as well as having a 3% case fatality rate. Those who survive the infection are thought to have long-term immunity.
- In the source country, the peak prevalence (i.e. maximum number of symptomatic people during the epidemic) was observed to be 8% of the population.
- Assume that the source country and your own country are epidemiologically equivalent (i.e. that disease transmission follows the same parameters in both settings).
Details rundown:
incubation period lasts 4 days on average = delta is 4 **-1
mean duration of symptoms before recovery or death is 5 days = gamma is 5 **-1
Symptomatic individuals are infectious, as well as having a 3% case fatality rate = cfr is .03
peak prevalence 8% of the population
mortality rate mu is not provided, but we can derive this from gamma and cfr; use this to indicate excess mortality for infected
infection rate beta is not provided, so we need to solve for beta
vaccine with 80% effectiveness available before epidemic begins
what is the minimum coverage for this vaccine?
Case fatality rate is mu divided by mu + gamma, so let’s calculate mu based on a CFR of 3% and a gamma of .20:
cfr = mu/(mu + gamma);
(cfr * mu) + (cfr * gamma) = mu;
(cfr * gamma) = mu - (cfr * mu);
(cfr * gamma) = mu * (1 - cfr);
(cfr * gamma)/(1 - cfr) = mu;
(.03 * .20)/(1 - .03) = mu;
(.006 / .97) = mu;
0.006185567 = mu;# 1 of 3
nicesubtitle <- "SEIR Model v1a1, respiratory virus with no vaccine c2w4_3"
print(nicesubtitle)## [1] "SEIR Model v1a1, respiratory virus with no vaccine c2w4_3"print("initial state values and parameters")## [1] "initial state values and parameters"# MODEL INPUTS:
N <- 100 # population size
duration <- 80 # total number of days
tsteps <- 0.1 # chunk in days
beta <- 4.1 # not provided -- infection rate day^-1
cfr <- .03 # case fatality rate
gamma <- 5 **-1 * (1 - cfr) # recovery rate adjusted to 3% cfr
mu <- 5 **-1 * (cfr) # excess mortality adjusted to 3% cfr
mu2 <- (cfr * gamma)/(1 - cfr) # excess mortality calculated as per above
delta <- 4 **-1 # incubation to infectious rate day^-1
R0 <- beta / gamma
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate
delta = delta, # incubation rate
mu = mu, # mortality rate, simple calc
mu2 = mu2, # mortality rate, derived
R0 = R0
)) ## beta gamma delta mu mu2 R0
## 4.10000 0.19400 0.25000 0.00600 0.00600 21.13402# initial_state_values and times
initial_state_values <- c(S = (N-1),
E = 0,
I = 1,
R = 0)
# TIMESTEPS:
times <- seq(from = 0, to = duration, by = tsteps)
# SEIR MODEL FUNCTION
seir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + E + I + R
lambda <- beta * I/N
dS <- -(lambda * S)
dE <- (lambda * S) -(delta * E)
dI <- (delta * E) -(gamma * I) -(mu * I)
dR <- (gamma * I)
return(list(c(dS, dE, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
output3 <- as.data.frame(ode(y = initial_state_values,
times = times,
func = seir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + E + I + R), digits=5) *100,
propor_Exp = round(E/(S + E + I + R), digits=5) *100,
preval_Inf = round(I/(S + E + I + R), digits=5) *100,
propor_Re = round(R/(S + E + I + R), digits=5) *100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + E + I + R)),
combo_E_I = E + I,
perc_E_I = round(combo_E_I/(S + E + I + R), digits=5) *100,
M_implied = N - (S + E + I + R)
)
print("peak Infected day, calculated: ")## [1] "peak Infected day, calculated: "max(output3$I)/N## [1] 0.3792863print("peak Infected day when I is at its max: ")## [1] "peak Infected day when I is at its max: "output3 %>%
arrange(-I, time) %>%
head(1)## time S E I R still_Su propor_Exp preval_Inf
## 1 9.3 0.1536014 30.51291 37.92863 30.46271 0.155 30.803 38.289
## propor_Re Reff combo_E_I perc_E_I M_implied
## 1 30.752 0.03277089 68.44154 69.092 0.9421456print("peak Infected day when preval_Inf % is at its max: ")## [1] "peak Infected day when preval_Inf % is at its max: "output3 %>%
arrange(-preval_Inf, time) %>%
head(1)## time S E I R still_Su propor_Exp preval_Inf
## 1 9.4 0.1312837 29.78158 37.92373 31.19851 0.133 30.072 38.293
## propor_Re Reff combo_E_I perc_E_I M_implied
## 1 31.502 0.02801584 67.70531 68.365 0.9649023print("peak Exposed & Infected day when (E + I) is at its max: ")## [1] "peak Exposed & Infected day when (E + I) is at its max: "output3 %>%
arrange(-(E+I), time) %>%
head(1)## time S E I R still_Su propor_Exp preval_Inf
## 1 7 4.635936 48.93116 31.53323 14.45269 4.657 49.151 31.675
## propor_Re Reff combo_E_I perc_E_I M_implied
## 1 14.518 0.9841588 80.46438 80.826 0.4469904print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output3 %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S E I R still_Su propor_Exp preval_Inf
## 1 7 4.635936 48.93116 31.53323 14.45269 4.657 49.151 31.675
## propor_Re Reff combo_E_I perc_E_I M_implied
## 1 14.518 0.9841588 80.46438 80.826 0.4469904print("last record for the run: ")## [1] "last record for the run: "output3 %>%
arrange(time) %>%
tail(1)## time S E I R still_Su propor_Exp
## 801 80 9.043305e-08 6.470519e-07 0.0001351099 96.99987 0 0
## preval_Inf propor_Re Reff combo_E_I perc_E_I M_implied
## 801 0 100 1.970324e-08 0.000135757 0 2.999996print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output3 %>%
select(-still_Su, -propor_Exp, -preval_Inf, -propor_Re, -Reff,
-combo_E_I, -perc_E_I, -M_implied) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line(linetype="dotdash") +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","chocolate","red","green")) +
scale_shape_manual(values = c(0,10,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", parameters['beta'],
" gamma of", round(parameters['gamma'],3),
" delta of", round(parameters['delta'],3),
" mu of", round(parameters['mu'],3),
"\n R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
- Create functions to seek appropriate beta for a peak infection rate of 8%
# 2 of 3
print("create a peak infection with closest value abs diff function, with some test runs")## [1] "create a peak infection with closest value abs diff function, with some test runs"# requires SEIR function from previous code chunk & prior parameters
# PEAKI FUNCTION (assumes seir_model model already loaded)
# parameters - vector with named elements for seir_model
SEIR_PEAKI <- function(parameters, peakval) {
# MODEL OUTPUT (solving the differential equations):
oderesult2 <- as.data.frame(ode(y = initial_state_values,
times = times,
func = seir_model,
parms = parameters))
peaki <- max(oderesult2$I)/N
# return absolute difference between calculated peak value and provided constant peakval
return(abs(peaki - peakval))
}
# test 1 with peak in a variable 'parm' and no peakval ref
SEIR_PEAKI(parameters =c(beta=4.1),peakval=0)## [1] 0.3792863# test 2 with peak in a variable 'parm' and no peakval ref
SEIR_PEAKI(parameters =c(beta=4.5),peakval=0)## [1] 0.3825844# test 3 with beta in a variable 'parm' and 8% peakval ref
SEIR_PEAKI(parameters =c(beta=4.1),peakval=0.08)## [1] 0.2992863# test 4 with beta in a variable 'parm' and 8% peakval ref
SEIR_PEAKI(parameters =c(beta=1.8),peakval=0.08)## [1] 0.2467476# Using optim() to find beta
# select upper limit of beta
betamax <- 80
# choose values to start your optimisation
beta_start <- .3
peakval <- 0.08
print("Using optim run 1 - No error code using L-BFGS-B method:")## [1] "Using optim run 1 - No error code using L-BFGS-B method:"# run optim ver1
(optimised1 <- optim(par = c(beta = beta_start)
, fn = SEIR_PEAKI
, peakval=peakval
, method = "L-BFGS-B"
))## $par
## beta
## 0.3826275
##
## $value
## [1] 3.682028e-08
##
## $counts
## function gradient
## 56 56
##
## $convergence
## [1] 0
##
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"print("Using optim run 2 - Warning on unreliability using Nelder-Mead method:")## [1] "Using optim run 2 - Warning on unreliability using Nelder-Mead method:"# run optim ver2
(optimised2 <- optim(par = c(beta = beta_start)
, fn = SEIR_PEAKI
, peakval=peakval
, method = "Nelder-Mead"
))## Warning in optim(par = c(beta = beta_start), fn = SEIR_PEAKI, peakval = peakval, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly## $par
## beta
## 0.3826275
##
## $value
## [1] 3.288581e-09
##
## $counts
## function gradient
## 48 NA
##
## $convergence
## [1] 0
##
## $message
## NULLprint("Using optimize run 3 - even with simplified function, does not appear to work:")## [1] "Using optimize run 3 - even with simplified function, does not appear to work:"# Using optimize() to find beta
# create hardocde 8 pct for single param optimize
SEIR_PEAK8 <- function(parameters) {
# MODEL OUTPUT (solving the differential equations):
oderesult2 <- as.data.frame(ode(y = initial_state_values,
times = times,
func = seir_model,
parms = parameters))
peaki <- max(oderesult2$I)/N
# return absolute difference between calculated peak value and hardcoded 8% constant peakval
return(abs(peaki - .08))
}
(optimised3 <- optimize(SEIR_PEAK8, c(0.01,betamax), tol = 0.0001))## $minimum
## [1] 79.99996
##
## $objective
## [1] 0.2992863# check results
# test 5 using optim run 1 - no error code using L-BFGS-B
print(paste0("test 5 using optim run 1, minimize diff from 8% peak value with beta = ", round(as.numeric(optimised1$par),6)))## [1] "test 5 using optim run 1, minimize diff from 8% peak value with beta = 0.382628"SEIR_PEAKI(parameters =c(beta= as.numeric(optimised1$par)), peakval=0.08)## [1] 3.682028e-08# test 6 using optim run 2 - warning using Nelder-Mead
print(paste0("test 6 using optim run 2, minimize diff from 8% peak value with beta = ", round(as.numeric(optimised2$par),6)))## [1] "test 6 using optim run 2, minimize diff from 8% peak value with beta = 0.382627"SEIR_PEAKI(parameters =c(beta= as.numeric(optimised2$par)), peakval=0.08)## [1] 3.288581e-09# test 7 using optimize run 3 - does not work?
print(paste0("test 7 using optimize run 3, minimize diff from 8% peak value with beta = ", round(as.numeric(optimised3$minimum),6)))## [1] "test 7 using optimize run 3, minimize diff from 8% peak value with beta = 79.999965"SEIR_PEAKI(parameters =c(beta= optimised3$minimum), peakval=0.08)## [1] 0.3285008# test best value and plot
nicesubtitle <- "SEIR Model v1a2, respiratory virus with no vaccine c2w4_3"
print(nicesubtitle)## [1] "SEIR Model v1a2, respiratory virus with no vaccine c2w4_3"(parameters <- c(
# beta = .32, #new beta
beta = as.numeric(optimised1$par), #new beta
gamma = gamma, # recovery rate
delta = delta, # incubation rate
mu = mu, # mortality rate, calculated
R0 = as.numeric(optimised1$par) / gamma # recalc with new beta
)) ## beta gamma delta mu R0
## 0.3826275 0.1940000 0.2500000 0.0060000 1.9723068# SEIR MODEL FUNCTION
seir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + E + I + R
lambda <- beta * I/N
dS <- -(lambda * S)
dE <- (lambda * S) -(delta * E)
dI <- (delta * E) -(gamma * I) -(mu * I)
dR <- (gamma * I)
return(list(c(dS, dE, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
output3 <- as.data.frame(ode(y = initial_state_values,
times = times,
func = seir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + E + I + R), digits=5) *100,
propor_Exp = round(E/(S + E + I + R), digits=5) *100,
preval_Inf = round(I/(S + E + I + R), digits=5) *100,
propor_Re = round(R/(S + E + I + R), digits=5) *100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + E + I + R)),
combo_E_I = E + I,
perc_E_I = round(combo_E_I/(S + E + I + R), digits=5) *100,
M_implied = N - (S + E + I + R)
)
print("peak Infected day, calculated: ")## [1] "peak Infected day, calculated: "max(output3$I)/N## [1] 0.08000004print("peak Infected day when I is at its max: ")## [1] "peak Infected day when I is at its max: "output3 %>%
arrange(-I, time) %>%
head(1)## time S E I R still_Su propor_Exp preval_Inf
## 1 47.8 48.42222 6.405023 8.000004 36.05757 48.968 6.477 8.09
## propor_Re Reff combo_E_I perc_E_I M_implied
## 1 36.464 0.9658052 14.40503 14.567 1.115183print("peak Infected day when preval_Inf % is at its max: ")## [1] "peak Infected day when preval_Inf % is at its max: "output3 %>%
arrange(-preval_Inf, time) %>%
head(1)## time S E I R still_Su propor_Exp preval_Inf propor_Re
## 1 47.9 48.27255 6.394692 8 36.21278 48.819 6.467 8.091 36.623
## Reff combo_E_I perc_E_I M_implied
## 1 0.9628668 14.39469 14.558 1.119983print("peak Exposed & Infected day when (E + I) is at its max: ")## [1] "peak Exposed & Infected day when (E + I) is at its max: "output3 %>%
arrange(-(E+I), time) %>%
head(1)## time S E I R still_Su propor_Exp preval_Inf
## 1 45.7 51.6643 6.574046 7.940363 32.80666 52.194 6.641 8.022
## propor_Re Reff combo_E_I perc_E_I M_implied
## 1 33.143 1.029423 14.51441 14.663 1.014639print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output3 %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S E I R still_Su propor_Exp preval_Inf
## 1 46.7 50.09772 6.505256 7.98297 34.35163 50.636 6.575 8.069
## propor_Re Reff combo_E_I perc_E_I M_implied
## 1 34.721 0.998691 14.48823 14.644 1.062422print("last record for the run: ")## [1] "last record for the run: "output3 %>%
arrange(time) %>%
tail(1)## time S E I R still_Su propor_Exp preval_Inf
## 801 80 24.61911 1.200934 2.228879 69.79254 25.162 1.227 2.278
## propor_Re Reff combo_E_I perc_E_I M_implied
## 801 71.332 0.4962768 3.429813 3.505 2.158532print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output3 %>%
select(-still_Su, -propor_Exp, -preval_Inf, -propor_Re, -Reff,
-combo_E_I, -perc_E_I, -M_implied) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line(linetype="dotdash") +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","chocolate","red","green")) +
scale_shape_manual(values = c(0,10,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", round(parameters['beta'],3),
" gamma of", round(parameters['gamma'],3),
" delta of", round(parameters['delta'],3),
" mu of", round(parameters['mu'],3),
"\n R0 of", round(parameters['R0'],3)),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")
- Calculate critical vaccine threshold for a perfectly effective vaccine:
pvacc_thresh = 1 - (gamma/beta);
pvacc_thresh = 1 - (.194/.3826275);
pvacc_thresh = 1 - .5070206;
pvacc_thresh = .4929794;
- Revise critical vaccine threshold for an 80% effective vaccine:
revised_pvacc = (1 - (gamma/beta))/veff;
revised_pvacc = (1 - (.194/.3826275))/0.80;
revised_pvacc = .6162244;
# 3 of 3
nicesubtitle <- "SEIR Model v1a3, respiratory virus with partially effective vaccine c2w4_3"
print(nicesubtitle)## [1] "SEIR Model v1a3, respiratory virus with partially effective vaccine c2w4_3"# revise parameters
beta <- as.numeric(parameters["beta"]) # best attempt
print("Critical Vaccine Threshold for perfect vaccine")## [1] "Critical Vaccine Threshold for perfect vaccine"(pvacc_thresh <- 1 - as.numeric((gamma/parameters["beta"])))## [1] 0.4929795R0 <- beta / gamma
veff <- 0.80 # 80% vaccine effectiveness
print("revised threshold for 80% effective vaccine")## [1] "revised threshold for 80% effective vaccine"(pvacc_thresh <- (1 - (1/R0))/veff) # recalc results to consider efficacy## [1] 0.6162244pvacc <- pvacc_thresh # reset to effective
(parameters <- c(
beta = beta, # infection rate
gamma = gamma, # recovery rate
delta = delta, # incubation rate
mu = mu, # mortality rate, simple calc
R0 = R0,
veff = veff,
pvacc = pvacc_thresh # reset to effective
)) ## beta gamma delta mu R0 veff pvacc
## 0.3826275 0.1940000 0.2500000 0.0060000 1.9723068 0.8000000 0.6162244# revise initial state values
# initial_state_values and times
initial_state_values <-
c(S = round((N-1) * (1-pvacc)),
E = 0,
I = 1,
R = round((N-1) * (pvacc)))
# SEIR MODEL FUNCTION
seir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
N <- S + E + I + R
lambda <- beta * I/N
dS <- -(lambda * S)
dE <- (lambda * S) -(delta * E)
dI <- (delta * E) -(gamma * I) -(mu * I)
dR <- (gamma * I)
return(list(c(dS, dE, dI, dR)))
})
}
# MODEL OUTPUT (solving the differential equations):
output3 <- as.data.frame(ode(y = initial_state_values,
times = times,
func = seir_model,
parms = parameters)) %>%
mutate(still_Su = round(S/(S + E + I + R), digits=5) *100,
propor_Exp = round(E/(S + E + I + R), digits=5) *100,
preval_Inf = round(I/(S + E + I + R), digits=5) *100,
propor_Re = round(R/(S + E + I + R), digits=5) *100,
Reff = (parameters['beta']/parameters['gamma'])
* (S/(S + E + I + R)),
combo_E_I = E + I,
perc_E_I = round(combo_E_I/(S + E + I + R), digits=5) *100,
M_implied = N - (S + E + I + R)
)
print("peak Infected day, calculated: ")## [1] "peak Infected day, calculated: "max(output3$I)/N## [1] 0.01print("peak Infected day when I is at its max: ")## [1] "peak Infected day when I is at its max: "output3 %>%
arrange(-I, time) %>%
head(1)## time S E I R still_Su propor_Exp preval_Inf propor_Re Reff combo_E_I
## 1 0 38 0 1 61 38 0 1 61 0.7494766 1
## perc_E_I M_implied
## 1 1 0print("peak Infected day when preval_Inf % is at its max: ")## [1] "peak Infected day when preval_Inf % is at its max: "output3 %>%
arrange(-preval_Inf, time) %>%
head(1)## time S E I R still_Su propor_Exp preval_Inf propor_Re Reff combo_E_I
## 1 0 38 0 1 61 38 0 1 61 0.7494766 1
## perc_E_I M_implied
## 1 1 0print("peak Exposed & Infected day when (E + I) is at its max: ")## [1] "peak Exposed & Infected day when (E + I) is at its max: "output3 %>%
arrange(-(E+I), time) %>%
head(1)## time S E I R still_Su propor_Exp preval_Inf propor_Re Reff combo_E_I
## 1 0 38 0 1 61 38 0 1 61 0.7494766 1
## perc_E_I M_implied
## 1 1 0print("point when R eff goes below 1: ")## [1] "point when R eff goes below 1: "output3 %>%
filter(Reff < 1) %>%
arrange(time) %>%
head(1)## time S E I R still_Su propor_Exp preval_Inf propor_Re Reff combo_E_I
## 1 0 38 0 1 61 38 0 1 61 0.7494766 1
## perc_E_I M_implied
## 1 1 0print("last record for the run: ")## [1] "last record for the run: "output3 %>%
arrange(time) %>%
tail(1)## time S E I R still_Su propor_Exp preval_Inf
## 801 80 35.72455 0.01941216 0.03001969 64.12924 35.759 0.019 0.03
## propor_Re Reff combo_E_I perc_E_I M_implied
## 801 64.191 0.7052803 0.04943185 0.049 0.09678061print("Plotting the proportion of people in each compartment over time")## [1] "Plotting the proportion of people in each compartment over time"output3 %>%
select(-still_Su, -propor_Exp, -preval_Inf, -propor_Re, -Reff,
-combo_E_I, -perc_E_I, -M_implied) %>%
melt(id = "time") %>%
mutate(proportion = value / sum(initial_state_values)) %>%
ggplot(aes(x=time, y=proportion, color=variable, shape=variable)) +
geom_line(linetype="dotdash") +
geom_jitter(show.legend = FALSE) +
scale_color_manual(values = c("blue","chocolate","red","green")) +
scale_shape_manual(values = c(0,10,4,1)) +
xlab("Time (days)") +
ylab("Proportion of people") +
labs(title=paste(" beta of", round(parameters['beta'],3),
" gamma of", round(parameters['gamma'],3),
" delta of", round(parameters['delta'],3),
" mu of", round(parameters['mu'],3),
"\n R0 of", round(parameters['R0'],3),
" veff of", round(parameters['veff'],3),
" pvacc of", round(parameters['pvacc'],3)
),
color = "Compartment",
subtitle= nicesubtitle) +
theme(legend.position="bottom")