A classic epidemiological model: the SIR model

Before diving into the real-life application, we first introduce the model that will be used.

There are many epidemiological models but we will use one of the most common one, the SIR model. The SIR model can be complexified to incorporate more specificities of the virus outbreak, but in this article we keep its simplest version. Tim Churches’ explanation of this model and how to fit it using R is so nice, I will reproduce it here with a few minor changes.

The basic idea behind the SIR model (Susceptible - Infectious - Recovered) of communicable disease outbreaks is that there are three groups (also called compartments) of individuals:

• S: those who are healthy but susceptible to the disease (i.e., at risk of being contaminated). At the start of the pandemic, S is the entire population since no one is immune to the virus.
• I: the infectious (and thus, infected) people
• R: individuals who were contaminated but who have either recovered or died. They are not infectious anymore.

These groups evolve over time as the virus progresses in the population:

• S decreases when individuals are contaminated and move to the infectious group I
• As people recover or die, they go from the infected group I to the recovered group R

To model the dynamics of the outbreak we need three differential equations to describe the rates of change in each group, parameterised by:

• \(\beta\), the infection rate, which controls the transition between S and I
• \(\gamma\), the removal or recovery rate, which controls the transition between I and R

Formally, this gives:

\[\frac{dS}{dt} = - \frac{\beta IS}{N} \text{ (Eq. 1)}\]

\[\frac{dI}{dt} = \frac{\beta IS}{N} - \gamma I \text{ (Eq. 2)}\]

\[\frac{dR}{dt} = \gamma I \text{ (Eq. 3)}\]

The first equation (Eq. 1) states that the number of susceptible individuals (S) decreases with the number of newly infected individuals, where new infected cases are the result of the infection rate (\(\beta\)) multiplied by the number of susceptible individuals (S) who had a contact with infectious individuals (I).

The second equation (Eq. 2) states that the number of infectious individuals (I) increases with the newly infected individuals (\(\beta I S\)), minus the previously infected people who recovered (i.e., \(\gamma I\) which is the removal rate \(\gamma\) multiplied by the infectious individuals I).

Finally, the last equation (Eq. 3) states that the recovered group (R) increases with the number of individuals who were infectious and who either recovered or died (\(\gamma I\)).

An epidemic develops as follows:

1. Before the start of the disease outbreak, S equals the entire population as no one has anti-bodies.
2. At the beginning of the outbreak, as soon as the first individual is infected, S decreases by 1 and I increases by 1 as well.
3. This first infectious individual contaminates (before recovering or dying) other individuals who were susceptible.
4. The dynamic continues, with recently contaminated individuals who in turn infect other susceptible people before they recover.

Visually, we have:

SIR model. Source: Kai Sasaki.

Before fitting the SIR model to the data, the first step is to express these differential equations as an R function, with respect to time t.

SIR <- function(time, state, parameters) {
  par <- as.list(c(state, parameters))
  with(par, {
    dS <- -beta * I * S / N
    dI <- beta * I * S / N - gamma * I
    dR <- gamma * I
    list(c(dS, dI, dR))
  })
}

Fitting a SIR model to the Belgium data

To fit the model to the data we need two things:

1. a solver for these differential equations
2. an optimiser to find the optimal values for our two unknown parameters, \(\beta\) and \(\gamma\)

The function ode() (for ordinary differential equations) from the {deSolve} R package makes solving the system of equations easy, and to find the optimal values for the parameters we wish to estimate, we can just use the optim() function built into base R.

Specifically, what we need to do is minimise the sum of the squared differences between \(I(t)\), which is the number of people in the infectious compartment \(I\) at time \(t\), and the corresponding number of cases as predicted by our model \(\hat{I}(t)\). This quantity is known as the residual sum of squares (RSS):

\[RSS(\beta, \gamma) = \sum_t \big(I(t) - \hat{I}(t) \big)^2\]

In order to fit a model to the incidence data for Belgium, we need a value N for the initial uninfected population. The population of Belgium in November 2019 was 11,515,793 people, according to Wikipedia. We will thus use N = 11515793 as the initial uninfected population.

Next, we need to create a vector with the daily cumulative incidence for Belgium, from February 4 (when our daily incidence data starts), through to March 30 (last available date at the time of publication of this article). We will then compare the predicted incidence from the SIR model fitted to these data with the actual incidence since February 4. We also need to initialise the values for N, S, I and R. Note that the daily cumulative incidence for Belgium is extracted from the {coronavirus} R package developed by Rami Krispin.

# devtools::install_github("RamiKrispin/coronavirus")
library(coronavirus)
data(coronavirus)

`%>%` <- magrittr::`%>%`

# extract the cumulative incidence
df <- coronavirus %>%
  dplyr::filter(country == "Belgium") %>%
  dplyr::group_by(date, type) %>%
  dplyr::summarise(total = sum(cases, na.rm = TRUE)) %>%
  tidyr::pivot_wider(
    names_from = type,
    values_from = total
  ) %>%
  dplyr::arrange(date) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(active = confirmed - death - recovered) %>%
  dplyr::mutate(
    confirmed_cum = cumsum(confirmed),
    death_cum = cumsum(death),
    recovered_cum = cumsum(recovered),
    active_cum = cumsum(active)
  )

# put the daily cumulative incidence numbers for Belgium from
# Feb 4 to March 30 into a vector called Infected
library(lubridate)

sir_start_date <- "2020-02-04"
sir_end_date <- "2020-03-30"

Infected <- subset(df, date >= ymd(sir_start_date) & date <= ymd(sir_end_date))$active_cum

# Create an incrementing Day vector the same length as our
# cases vector
Day <- 1:(length(Infected))

# now specify initial values for N, S, I and R
N <- 11515793
init <- c(
  S = N - Infected[1],
  I = Infected[1],
  R = 0
)

Note that what is needed are currently infected persons (cumulative infected minus the removed, i.e. recovered or dead). However, numbers of recovered persons are hard to obtain and probably underestimated due to underreporting bias. I thus consider the cumulative number of infected people, which is probably not an issue here since the number of recovered cases is negligible at the time of the analysis.

Then we need to define a function to calculate the RSS, given a set of values for \(\beta\) and \(\gamma\).

# define a function to calculate the residual sum of squares
# (RSS), passing in parameters beta and gamma that are to be
# optimised for the best fit to the incidence data
RSS <- function(parameters) {
  names(parameters) <- c("beta", "gamma")
  out <- ode(y = init, times = Day, func = SIR, parms = parameters)
  fit <- out[, 3]
  sum((Infected - fit)^2)
}

Finally, we can fit the SIR model to our data by finding the values for \(\beta\) and \(\gamma\) that minimise the residual sum of squares between the observed cumulative incidence (observed in Belgium) and the predicted cumulative incidence (predicted by our model). We also need to check that our model has converged, as indicated by the message shown below:

# now find the values of beta and gamma that give the
# smallest RSS, which represents the best fit to the data.
# Start with values of 0.5 for each, and constrain them to
# the interval 0 to 1.0

# install.packages("deSolve")
library(deSolve)

Opt <- optim(c(0.5, 0.5),
  RSS,
  method = "L-BFGS-B",
  lower = c(0, 0),
  upper = c(1, 1)
)

# check for convergence
Opt$message

## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

Convergence is confirmed. Note that you may find different estimates for different choices of initial values or constraints. This proves that the fitting process is not stable. Here is a potential solution for a better fitting process.

Now we can examine the fitted values for \(\beta\) and \(\gamma\):

Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par

##      beta     gamma 
## 0.5841185 0.4158816

Remember that \(\beta\) controls the transition between S and I (i.e., susceptible and infectious) and \(\gamma\) controls the transition between I and R (i.e., infectious and recovered). However, those values do not mean a lot but we use them to get the fitted numbers of people in each compartment of our SIR model for the dates up to March 30 that were used to fit the model, and compare those fitted values with the observed (real) data.

# time in days for predictions
t <- 1:as.integer(ymd(sir_end_date) + 1 - ymd(sir_start_date))

# get the fitted values from our SIR model
fitted_cumulative_incidence <- data.frame(ode(
  y = init, times = t,
  func = SIR, parms = Opt_par
))

# add a Date column and the observed incidence data
library(dplyr)
fitted_cumulative_incidence <- fitted_cumulative_incidence %>%
  mutate(
    Date = ymd(sir_start_date) + days(t - 1),
    Country = "Belgium",
    cumulative_incident_cases = Infected
  )

# plot the data
library(ggplot2)
fitted_cumulative_incidence %>%
  ggplot(aes(x = Date)) +
  geom_line(aes(y = I), colour = "red") +
  geom_point(aes(y = cumulative_incident_cases), colour = "blue") +
  labs(
    y = "Cumulative incidence",
    title = "COVID-19 fitted vs observed cumulative incidence, Belgium",
    subtitle = "(Red = fitted from SIR model, blue = observed)"
  ) +
  theme_minimal()

From the above graph we see that the number of observed confirmed cases follows (unfortunately) the number of confirmed cases expected by our model. The fact that both trends are overlapping indicates that the pandemic is clearly in an exponential phase in Belgium. More data would be needed to see whether this trend is confirmed in the long term.

The following graph is similar than the previous one, except that the y-axis is measured on a log scale. This kind of plot is called a semi-log plot or more precisely a log-linear plot because only the y-axis is transformed with a logarithm scale. Transforming the scale in log has the advantage that it is more easily readable in terms of difference between the observed and expected number of confirmed cases and it also shows how the number of observed confirmed cases differs from an exponential trend.

fitted_cumulative_incidence %>%
  ggplot(aes(x = Date)) +
  geom_line(aes(y = I), colour = "red") +
  geom_point(aes(y = cumulative_incident_cases), colour = "blue") +
  labs(
    y = "Cumulative incidence",
    title = "COVID-19 fitted vs observed cumulative incidence, Belgium",
    subtitle = "(Red = fitted from SIR model, blue = observed)"
  ) +
  theme_minimal() +
  scale_y_log10(labels = scales::comma)

The plot indicates that, at the beginning of the pandemic and until March 12, the number of confirmed cases stayed below what would be expected in an exponential phase. In particular, the number of confirmed cases stayed constant at 1 case from February 4 to February 29. From March 13 and until March 30, the number of confirmed cases kept increasing at a rate close to an exponential rate.

We also notice a small jump between March 12 and March 13, which may potentially indicate an error in the data collection, or a change in the testing/screening methods.

Reproduction number \(R_0\)

Our SIR model looks like a good fit to the observed cumulative incidence data in Belgium, so we can now use our fitted model to calculate the basic reproduction number \(R_0\), also referred as basic reproduction ratio, and which is closely linked to \(\beta\) and \(\gamma\).²

The basic reproduction number \(R_0\) gives the average number of susceptible people who are infected by each infectious person where all individuals are susceptible to infection. In other words, the reproduction number refers to the number of healthy people that get infected per number of sick people. When \(R_0 > 1\) the disease starts spreading in a population, but not if \(R_0 < 1\). Usually, the larger the value of \(R_0\), the harder it is to control the epidemic and the higher the probability of a pandemic.

Formally, we have:

\[R_0 = \frac{\beta}{\gamma}\]

We can compute it in R:

Opt_par

##      beta     gamma 
## 0.5841185 0.4158816

R0 <- as.numeric(Opt_par[1] / Opt_par[2])
R0

## [1] 1.404531

An \(R_0\) of 1.4 is below values found by others for COVID-19 and the \(R_0\) for SARS and MERS, which are similar diseases also caused by coronavirus. Furthermore, in the literature, it has been estimated that the reproduction number for COVID-19 is approximately 2.7 (with \(\beta\) close to 0.54 and \(\gamma\) close to 0.2). Our reproduction number being lower is mainly due to the fact that the number of confirmed cases stayed constant and equal to 1 at the beginning of the pandemic.

A \(R_0\) of 1.4 means that, on average in Belgium, 1.4 persons are infected for each infected person.

For simple models, the proportion of the population that needs to be effectively immunized to prevent sustained spread of the disease, known as the “herd immunity threshold”, has to be larger than \(1 - \frac{1}{R_0}\) (Fine, Eames, and Heymann 2011).

The reproduction number of 1.4 we just calculated suggests that, given the formula 1 - (1 / 1.4), 28.8% of the population should be immunized to stop the spread of the infection. With a population in Belgium of approximately 11.5 million, this translates into roughly 3.3 million people.

Using our model to analyze the outbreak if there was no intervention

It is instructive to use our model fitted to the first 56 days of available data on confirmed cases in Belgium, to see what would happen if the outbreak were left to run its course, without public health intervention.

# time in days for predictions
t <- 1:120

# get the fitted values from our SIR model
fitted_cumulative_incidence <- data.frame(ode(
  y = init, times = t,
  func = SIR, parms = Opt_par
))

# add a Date column and join the observed incidence data
fitted_cumulative_incidence <- fitted_cumulative_incidence %>%
  mutate(
    Date = ymd(sir_start_date) + days(t - 1),
    Country = "Belgium",
    cumulative_incident_cases = c(Infected, rep(NA, length(t) - length(Infected)))
  )

# plot the data
fitted_cumulative_incidence %>%
  ggplot(aes(x = Date)) +
  geom_line(aes(y = I), colour = "red") +
  geom_line(aes(y = S), colour = "black") +
  geom_line(aes(y = R), colour = "green") +
  geom_point(aes(y = cumulative_incident_cases),
    colour = "blue"
  ) +
  scale_y_continuous(labels = scales::comma) +
  labs(y = "Persons", title = "COVID-19 fitted vs observed cumulative incidence, Belgium") +
  scale_colour_manual(name = "", values = c(
    red = "red", black = "black",
    green = "green", blue = "blue"
  ), labels = c(
    "Susceptible",
    "Recovered", "Observed", "Infectious"
  )) +
  theme_minimal()

The same graph in log scale for the y-axis and with a legend for better readability:

# plot the data
fitted_cumulative_incidence %>%
  ggplot(aes(x = Date)) +
  geom_line(aes(y = I, colour = "red")) +
  geom_line(aes(y = S, colour = "black")) +
  geom_line(aes(y = R, colour = "green")) +
  geom_point(aes(y = cumulative_incident_cases, colour = "blue")) +
  scale_y_log10(labels = scales::comma) +
  labs(
    y = "Persons",
    title = "COVID-19 fitted vs observed cumulative incidence, Belgium"
  ) +
  scale_colour_manual(
    name = "",
    values = c(red = "red", black = "black", green = "green", blue = "blue"),
    labels = c("Susceptible", "Observed", "Recovered", "Infectious")
  ) +
  theme_minimal()

fit <- fitted_cumulative_incidence

# peak of pandemic
fit[fit$I == max(fit$I), c("Date", "I")]

##          Date        I
## 89 2020-05-02 531000.4

# severe cases
max_infected <- max(fit$I)
max_infected * 0.2

## [1] 106200.1

# cases with need for intensive care
max_infected * 0.06

## [1] 31860.03

# deaths with supposed 4.5% fatality rate
max_infected * 0.045

## [1] 23895.02

Ascertainment rates

In the previous analyses and graphs, it is assumed that the number of confirmed cases represent all the cases that are infectious. This is far from reality as only a proportion of all cases are screened, detected and counted in the official figures. This proportion is known as the ascertainment rate.

The ascertainment rate is likely to vary during the course of an outbreak, in particular if testing and screening efforts are increased, or if detections methods are changed. Such changing ascertainment rates can be easily incorporated into the model by using a weighting function for the incidence cases.

In his first article, Tim Churches demonstrates that a fixed ascertainment rates of 20% makes little difference to the modelled outbreak with no intervention, except that it all happens a bit more quickly.

More sophisticated models

More sophisticated models could also be used to better reflect real-life transmission processes. For instance, another classical model in disease outbreak is the SEIR model. This extended model is similar to the SIR model, where S stands for Susceptible and R stands for Recovered, but the infected people are divided into two compartments:

1. E for the Exposed/infected but asymptomatic
2. I for the Infected and symptomatic

These models belong to the continuous-time dynamic models that assume fixed transition rates. There are other stochastic models that allow for varying transition rates depending on attributes of individuals, social networking, etc.

Modelling the epidemic trajectory using log-linear models

As noted above, the initial exponential phase of an outbreak, when shown in a log-linear plot (the y-axis on a log scale and the x-axis without transformation), appears (somewhat) linear. This suggests that we can model epidemic growth, and decay, using a simple log-linear model of the form:

\[log(y)=rt+b\]

where y is the incidence, r is the growth rate, t is the number of days since a specific point in time (typically the start of the outbreak), and b is the intercept. In this context, two log-linear models:

1. one to the growth phase (before the peak), and
2. one to the decay phase (after the peak)

are fitted to the epidemic (incidence cases) curve.

The doubling and halving time estimates which you very often hear in the news can be estimated from these log-linear models. Furthermore, these log-linear models can also be used on the epidemic trajectory to estimate the reproduction number \(R_0\) in the growth and decay phases of the epidemic.

The {incidence} package in R, part of the R Epidemics Consortium (RECON) suite of packages for epidemic modelling and control, makes the fitting of this kind of models very convenient.

Estimating changes in the effective reproduction number \(R_e\)

In our model, we set a reproduction number \(R_0\) and kept it constant. It would nonetheless be useful to estimate the current effective reproduction number \(R_e\) on a day-by-day basis so as to track the effectiveness of public health interventions, and possibly predict when an incidence curve will start to decrease.

The {EpiEstim} package in R can be used to estimate \(R_e\) and allow to take into consideration human travel from other geographical regions in addition to local transmission (Cori et al. 2013; Thompson et al. 2019).

More sophisticated projections

In addition to naïve predictions based on a simple SIR model, more advanced and complex projections are also possible, notably, with the {projections} package. This packages uses data on daily incidence, the serial interval and the reproduction number to simulate plausible epidemic trajectories and project future incidence.