basic-SIR.R

#Here is a basic reproduction of the SIR model with a level of infection rate of 0.5
#and a recovery rate of 0.25 replicated for 100 days. 


# LIBRARY--------------

library(deSolve)
library(reshape2)
library(ggplot2)

# BUILDIN UP A BASIC SIR MODEL


# INPUTS---------

N<- 100000                  # population

state_values<- c(S = N -1,  # susceptible
                I = 1,      # infected
                R = 0)      # recovered

parameters<- c(beta = 1/2,  # infection rate
               gamma = 1/4) # recovery rate

              
# TIMEFRAME---------

times<- seq(0, 100, by = 1)

# MODEL----------

sir_model<- function(time, state, parameters){
  with(as.list(c(state, parameters)),{
    N<- S + I + R
    lambda = beta * I/N    # force of infection
    dS<- - lambda * S 
    dI<- lambda * S - gamma * I
    dR<- gamma * I
    return(list(c(dS,dI,dR)))
  })
}

# OUTPUT---------

 # calculation of the differential equations

output<- as.data.frame(ode(y = state_values,
                           times = times,
                           func = sir_model,
                           parms = parameters))

output_full<- melt(as.data.frame(output), id = "time")

# the proportion of each compartment ( S, I, R)

output_full$proportion<- output_full$value/sum(state_values)
                              
                           
# PLOT---------

# plot 1: the prevalence of infection

ggplot(data = output, aes(x = time, y = I)) + 
  geom_line() + 
  xlab("Time(days)") +
  ylab("Number of Infected") + 
  labs("SIR Model: prevalence of infection")

# plot 2: the SIR model 

ggplot(output_full, aes(time, proportion, color = variable, group = variable)) + 
  geom_line() + 
  xlab("Time(days)") +
  ylab("Prevalence") + 
  labs(color = "Compartment", title = "SIR Model")

# Reff---------------------------

# calculation of the effective reproduction number (Reff)

output$reff<- parameters["beta"]/parameters["gamma"] *
  output$S/(output$S+output$I+output$R)

# calculation of the reproduction number (R0)

output$R0<- parameters["beta"]/parameters["gamma"]

# PLOT Reff------

# plot 3: levels of reproduction of infection

ggplot() + 
  geom_line (data = output, aes(x = time, y = reff)) + 
  geom_point (data = output, aes(x = time, y = reff), colour = "red") + 
  geom_line (data = output, aes(x = time, y = R0), colour = "green") + 
  xlab("Time(days)") + 
  ylab("Reff") + 
  labs(title = paste("Reproduction number levels with: Beta = ", 
             parameters["beta"], " and Gamma = ", parameters["gamma"]))