## 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.