Skeleton code for the DGHS Workshop

intro_sir.R

## Import packages
library(deSolve)

## Basic SIR using number of people (not proportion of pop)
SIR <- function(t, x, parms){
    with(as.list(c(parms, x)), {
        N = S + I + R
        dS <- -beta*S*I/N
        dI <- beta*S*I/N - r*I
        dR <- r*I
        result <- c(dS, dI, dR)
        list(result)
    }) 
}

## Define time and initial state conditions
inits <- c(S = 999999, I = 1, R = 0)   # closed population, 1m people
dt <- seq(0, 60, .5)                 # 60 days, hourly increments
beta <- .75 * (12/7)                   # given -- 12 contacts per week
r <- 1/7                               # duration = 1 week

SIR_sim <- as.data.frame(ode(inits, dt, SIR, parms=c(beta, r)))
matplot(SIR_sim$time, SIR_sim[, 2:4], type = "l", lty = 1)