## Introduction to Infectious Disease Dynamics (Chain Binomial - Tutorial 3)
## Clinic on the Meaningful Modeling of Epidemiological Data
## African Institute for Mathematical Sciences, Muizenberg, RSA
##
## Juliet R.C. Pulliam, 2012-2015
##
## Tutorial 3, Benchmark question 3 for Math Models in Med \& PH course
rm(list=ls()) # Clear all variables and functions
## SIR.CB() -- Function that runs a chain binomial SIR model for a given population
## size and R0 value, outputting a data frame with columns representing, time,
## the number of susceptibles at each time, and the number of cases at each time.
## The chain binomial is a stochastic model, so the output varies between function
## calls, even for the same parameter values and initial conditions.
sir.cb <- function(R0,N,MAXTIME,I0=1,S0=N-I0){
qq <- 1-R0/N # Pairwise probability of avoiding potentially infectious contact
Cases <- I0
Sus <- S0
for(Time in 1:MAXTIME){
Cases <- c(Cases,rbinom(1,Sus[Time],(1-qq^Cases[Time])))
Sus <- c(Sus,Sus[Time]-Cases[Time+1])
}
return(data.frame(Time=0:MAXTIME,Cases,Sus))
}
# Run the function once, for an R0 of 1.5, a population size of 1000, and 60 time
# steps
epi <- sir.cb(1.5,1000,60)
# Examine the resulting data frame
head(epi)
# Plot the number of cases through time and label the plot appropriately
par(mar=c(5,5,1,1))
plot(0:60,epi$Cases,
type="s", # Use a 'step' plot because time is treated as discrete
bty="n",
lwd=3,
cex.lab=2,
ylim=c(0,100),
cex.axis = 1.1,
xlab="Time",
ylab="Infected")
text(40,80,expression(R[0]==1.5),cex=2)
# PLOT.CB -- Function that runs SIR.CB for specified parameter values, plots
# the number of cases through time (if plot==T), and returns the vector of values
# of the number of cases through time
plot.cb <- function(R0,N,MAXTIME=60,lwd=1,col="grey",plot=T){
cases <- sir.cb(R0,N,MAXTIME)$Cases
if(plot) lines(0:MAXTIME,cases,type="s",lwd=lwd,col=col)
return(cases)
}
plot.cb(1.5,1000,60)
# Set up an empty plot with pre-labelled axes
par(mar=c(5,5,1,1))
plot(0:60,epi$Cases,
type="s", # Use a 'step' plot because time is treated as discrete
bty="n",
lwd=0,
cex.lab=2,
ylim=c(0,100),
cex.axis = 1.1,
xlab="Time",
ylab="Infected")
text(40,80,expression(R[0]==1.5),cex=2) # Add the R0 value used to the plot
# Call plot.cb() 3 times to plot 3 runs of the SIR chain binomial; do this
# using the built-in replicate() function, which automatically concatenates the
# output into a matrix
runs <- replicate(3,plot.cb(1.5,1000,60))
# Examine the upper left portion of the stored matrix: each column represents
# a different run, and each row represents a time point in the run
head(runs)
## Calculate the average number of cases in each timestep. Do this by using
## the apply() function to apply the function mean() to the rows of the matrix,
## runs, created above:
ave.case.t <- apply(runs,1,mean)
## Now, complete the following tasks on your own.
##
## Set up an empty plot with pre-labelled axes, just like before:
# Add the R0 value used to the plot:
## Call plot.cb() 300 times to plot 300 runs of the SIR chain binomial. Again,
## do this using the built-in replicate() function to automatically concatenate
## the output:
## Calculate the average number of cases in each timestep:
## Calculate the median number of cases in each timestep:
## Add lines to your plot that represent the mean and median number of cases
## through time. Use a thick red line for the mean values and a thick blue line
## for the median values:
## Save your plot as a PDF, and email it to Steve.