SIR model with a constant force of infection

library(deSolve)   # package to solve the model
library(reshape2)  # package to change the shape of the model output
library(ggplot2)   # package for plotting

#1.1 The input data from the instructions were as follows:

initial_state_values <- c(S = 999999,  # the whole population we are modelling is susceptible to infection
                          I = 1,       # the epidemic starts with a single infected person
                          R = 0)       # there is no prior immunity in the population

# Vector storing the parameters describing the transition rates in units of days^-1
parameters <- c(lambda = 0.2,  # the force of infection, which acts on susceptibles
                gamma = 0.071)   # the rate of recovery, which acts on those infected

1 TIMESTEPS:

2 Vector storing the sequence of timesteps to solve the model at

times <- seq(from = 0, to = 60, by = 1)   # from 0 to 60 days in daily intervals

# SIR MODEL FUNCTION: 

# The model function takes as input arguments (in the following order): time, state and parameters
sir_model <- function(time, state, parameters) {  

    with(as.list(c(state, parameters)), {  # tell R to unpack variable names from the state and parameters inputs
        
    # The differential equations
      dS <- -lambda * S               # people move out of (-) the S compartment at a rate lambda (force of infection)
      dI <- lambda * S - gamma * I    # people move into (+) the I compartment from S at a rate lambda, 
                                      # and move out of (-) the I compartment at a rate gamma (recovery)
      dR <- gamma * I                 # people move into (+) the R compartment from I at a rate gamma
      
    # Return the number of people in the S, I and R compartments at each timestep 
    # (in the same order as the input state variables)
    return(list(c(dS, dI, dR))) 
    })
  
}

3 MODEL OUTPUT (solving the differential equations):

4 Solving the differential equations using the ode integration algorithm

output <- as.data.frame(ode(y = initial_state_values, 
                            times = times, 
                            func = sir_model,
                            parms = parameters))
# Printing the model output returns a dataframe with columns time (containing the times vector), 
# S (containing the number of susceptible people at each timestep),
# I (containing the number of infected people at each timestep) and 
# R (containing the number of recovered people at each timestep).

#print result

output

##    time            S         I         R
## 1     0 9.999990e+05      1.00      0.00
## 2     1 8.187299e+05 174777.72   6492.35
## 3     2 6.703194e+05 305893.91  23786.71
## 4     3 5.488111e+05 402084.91  49104.01
## 5     4 4.493285e+05 470446.31  80225.18
## 6     5 3.678790e+05 516735.12 115385.84
## 7     6 3.011939e+05 545615.87 153190.26
## 8     7 2.465967e+05 560862.16 192541.17
## 9     8 2.018963e+05 565521.16 232582.57
## 10    9 1.652987e+05 562048.24 272653.08
## 11   10 1.353351e+05 552416.97 312247.95
## 12   11 1.108030e+05 538208.84 350988.19
## 13   12 9.071779e+04 520686.85 388595.36
## 14   13 7.427343e+04 500855.36 424871.21
## 15   14 6.080993e+04 479508.97 459681.09
## 16   15 4.978695e+04 457272.52 492940.52
## 17   16 4.076210e+04 434633.55 524604.34
## 18   17 3.337318e+04 411968.87 554657.94
## 19   18 2.732365e+04 389566.18 583110.17
## 20   19 2.237071e+04 367641.60 609987.69
## 21   20 1.831558e+04 346354.03 635330.39
## 22   21 1.499553e+04 325816.72 659187.75
## 23   22 1.227730e+04 306106.74 681615.96
## 24   23 1.005180e+04 287272.55 702675.65
## 25   24 8.229717e+03 269340.26 722430.02
## 26   25 6.737922e+03 252318.55 740943.52
## 27   26 5.516543e+03 236202.75 758280.70
## 28   27 4.516563e+03 220978.03 774505.41
## 29   28 3.697848e+03 206622.00 789680.15
## 30   29 3.027542e+03 193106.82 803865.64
## 31   30 2.478741e+03 180400.79 817120.47
## 32   31 2.029421e+03 168469.69 829500.89
## 33   32 1.661549e+03 157277.79 841060.66
## 34   33 1.360361e+03 146788.67 851850.97
## 35   34 1.113769e+03 136965.81 861920.42
## 36   35 9.118772e+02 127773.10 871315.03
## 37   36 7.465818e+02 119175.15 880078.27
## 38   37 6.112494e+02 111137.59 888251.16
## 39   38 5.004486e+02 103627.26 895872.29
## 40   39 4.097326e+02  96612.31 902977.95
## 41   40 3.354606e+02  90062.30 909602.24
## 42   41 2.746519e+02  83948.23 915777.12
## 43   42 2.248659e+02  78242.58 921532.55
## 44   43 1.841046e+02  72919.29 926896.61
## 45   44 1.507321e+02  67953.71 931895.56
## 46   45 1.234090e+02  63322.64 936553.95
## 47   46 1.010387e+02  59004.19 940894.77
## 48   47 8.272348e+01  54977.82 944939.46
## 49   48 6.772825e+01  51224.20 948708.07
## 50   49 5.545119e+01  47725.23 952219.32
## 51   50 4.539959e+01  44463.92 955490.68
## 52   51 3.717003e+01  41424.38 958538.45
## 53   52 3.043225e+01  38591.73 961377.84
## 54   53 2.491581e+01  35952.05 964023.04
## 55   54 2.039934e+01  33492.32 966487.29
## 56   55 1.670156e+01  31200.38 968782.92
## 57   56 1.367408e+01  29064.88 970921.44
## 58   57 1.119539e+01  27075.22 972913.58
## 59   58 9.166008e+00  25221.49 974769.34
## 60   59 7.504492e+00  23494.46 976498.03
## 61   60 6.144157e+00  21885.51 978108.35

#Plotting the output:

output_long <- melt(as.data.frame(output), id = "time")                  # turn output dataset into long format

ggplot(data = output_long,                                               # specify object containing data to plot
       aes(x = time, y = value, colour = variable, group = variable)) +  # assign columns to axes and groups
  geom_line() +                                                          # represent data as lines
  xlab("Time (days)")+                                                   # add label for x axis
  ylab("Number of people") +                                             # add label for y axis
  labs(colour = "Compartment")                                           # add legend title

1.1.1 Based on the plot, describe the pattern of the epidemic over the 2 month period. How does the number of people in the susceptible, infected and recovered compartment change over time? After how many days does the epidemic reach its peak? After how many days does it end?

The number of infected people quickly increases, reaching a peak of 500000 infected people after around 7 days, before steadily decreasing again. The number of recovered people starts to rise shortly after the first people become infected. It increases steadily (but less sharply than the curve of infected people) until the whole population has become immune - by day 53, 99% are in the R compartment, and nearly no susceptible people remain after 60 days.