SIR Model
Carlos C Vera Recio
12/30/2021
An example of a simple SIR model in R
- The SIR model is the epidemiology workhorse
- SIR refers to “Susceptible”, “Infected” and “Recovered”
- The SIR model can be represented using the following figure:
Figure 1
- We can model the change of each the “Susceptible”, “Infected” and “Recovered” compartments over time using the following differential equations:
Figure 2
- Here is a summary of the parameters in Figure 2
\(\lambda\) = force of infection. We assume it changes over time and that is depends on the fraction of the population that is infected.
\(\beta\) = rate of infection
\(\gamma\) = rate of recovery
- I use the deSolve package to solve this model in R
- We will also use the reshape2 and ggplot2 packages later on
- First, let’s load all packages
# load libraries
library(deSolve)
library(ggplot2)
library(reshape2)- We will now do an example of an SIR model
Our population will have 10^6 members, all susceptible to an incoming infection
One members of the population contracts the disease
The recovery time is, on average, 10 days
The infection rate is 1
We will monitor this population for 60 days
- Let’s prepare these inputs
- We will begin by creating a vector of values that represents the initial number of people in the Susceptible, Infected and Recovered compartments.
# initial population values
initial_number_susceptible = 10e6 - 1
initial_number_recovered = 0
initial_number_infected = 1
initial_state_values <- c(
S = initial_number_susceptible,
I = initial_number_infected,
R = initial_number_recovered
)- We will create a parameters vector that includes the values for \(\beta\) and \(\gamma\)
# model parameters
infection_rate = 1
recovery_rate = 0.1
parameters <- c(
beta = infection_rate,
gamma = recovery_rate
)- We will create a variable that represents how often and how long we will collect data on the infectious process.
- We will follow the data daily for 60 days
# time steps
times = seq(
from = 0,
to = 60,
by = 1
)- We will now write a function that represents the differential equations in figure 2
# SIR MODEL FUNCTION
sir_model <- function(
time,
state,
parameters
)
{
with(as.list(c(state, parameters)), {
# Differential equations
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))) # return output
})
}- Next we run the model and collect the output
# run the model
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters))- Finally, let’s plot the data!
# reformat data frame
output_long <- melt(as.data.frame(output), id = "time")
# Plot the number of people in each compartment over time
gp <- ggplot(data = output_long,
aes(x = time, y = value, colour = variable, group = variable)) +
geom_line() +
xlab("Time (days)")+
ylab("Number of people") +
labs(title = paste("Susceptible, Infected and Recovered"))
gp
That’s it! As you can see, this was actually a pretty simple exercise. You can adjust the initial values and model parameters so that the model applies for your particular case of interest (for example, the COVID-19 pandemic). This is a simple implementation of the SIR model that does not accomodate population turnover, effect of vaccination or waning immunity. As you might imagine, these assumptions do not apply in many scenarios. So let’s do it!. Next, we’ll code several more complex versions of the SIR model that can accommodate the aforementioned factors.
SIR model that accomodates population turnover
- The next SIR model we’ll discuss be represented using the following figure:
Figure 3
\(\lambda\) = force of infection. We assume it changes over time and that is depends on the fraction of the population that is infected.
\(\beta\) = rate of infection
\(\gamma\) = rate of recovery
\(\mu\) = rate of people exiting population (all cause mortality)
\(b\) = rate of people entering compartment (birth rate)
All right! Let’s code it up!
# initial population values
initial_number_susceptible = 10e6 - 1
initial_number_recovered = 0
initial_number_infected = 1
initial_population = initial_number_infected+initial_number_recovered+initial_number_susceptible
initial_state_values <- c(
S = initial_number_susceptible,
I = initial_number_infected,
R = initial_number_recovered,
N = initial_population
)
# time steps
times <- seq(from = 0, to = 100, by = 1/365) # from 0 to 100 YEARS in DAILY intervals
# model parameters
infection_rate = 1*365
recovery_rate = 0.1*365
death_rate = 1/(70)
birth_rate = death_rate
parameters <- c(
beta = infection_rate,
gamma = recovery_rate,
b = birth_rate,
miu = death_rate
)
# SIR MODEL FUNCTION
sir_model <- function(
time,
state,
parameters
)
{
with(as.list(c(state, parameters)), {
# Differential equations
N <- S+I+R
lambda <- beta*(I/N)
dS <- -lambda*S + b*N - miu*S
dI <- lambda*S - gamma*I - miu*I
dR <- gamma*I - R*miu
return(list(c(dS, dI, dR, N))) # return output
})
}
# run the model
output <- as.data.frame(ode(y = initial_state_values,
times = times,
func = sir_model,
parms = parameters))
# reformat data frame
plot.cols = c('time', 'S','I','R')
output_long <- melt(as.data.frame(output[,plot.cols]), id = "time")
# Plot the proportion of people in each compartment over time
gp1 <- ggplot(data = output_long,
aes(x = time, y = value, colour = variable, group = variable)) +
geom_line() +
xlab("Time (years)")+
ylab("Number of people") +
labs(title = paste("Infected over time"))
gp1
As a new surplus of susceptible individuals enter the population, oscillations in the number of infected people occur throughout history. This is represented in the model through the terms that we added to the model representing population turnover. For this example, we assumed birth and death are in equal rate. In reality populations tend to increase over time. This can, of course, be adjusted through the death and birth rate parameters (\(b\) and \(\mu\)).
One last important caveat is that the parameters in this model, for example \(\gamma\) and \(\beta\), need to be estimated from the population. Poor estimation of these values would lead to a poor performance by the SIR model. Eventually, we’ll include in this R-markdown how estimation of these parameters can be performed from a data set.