SIMC_SIR_Exercise_1.Rmd

---
title: "SIMC_SIR_Exercise_1"
author: "Daniel T Citron"
date: '2023-07-09'
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Import Libraries

```{r Download Libraries}
library(deSolve)   # numerical integration package
library(ggplot2)   # plotting package
library(tidyverse) # packages for handling and reshaping data
library(reshape2) 
library(magrittr)
library(stats) # statistics package
```

# Exercise 1: Numerical solution to the SIR model

We will use a numerical ODE solver to integrate the SIR model. Our state variables are (S, I, R) and our parameters are beta and gamma.

```{r Define SIR model function to be integrated}
sir_model <- function(time, state, parameters) {
  # Unpack Parameters
  with(as.list(c(state,parameters)), {
      # Calculate total population size
      N <-  S + I + R
      
      # Calculate force of infection
      lambda <- beta * I / N
      
      # Calculate derivatives
      dS <- -lambda * S
      dI <- lambda * S - gamma * I
      dR <- gamma * I
      
      # Return derivative
      return(list(c(dS, dI, dR)))
    }
  )
}
```

Define the initial state variables, parameters, and time span. Integrate the SIR function defined in the previous block.
```{r}
# Initial values for state variables
initial_state_values <- c(S = 9999,
                         I = 1,
                         R = 0
)

# Parameter values
parameters <-  c(beta = 1, # per-contact infection rate
                 gamma = .5 # recovery rate
)

# Define a list of time steps 
times = seq(from = 0, to = 50, by = .1)

# Integrate
sir.output <- as.data.frame(ode(y = initial_state_values, 
                            times = times, 
                            func = sir_model,
                            parms = parameters))
```

Plot the output:
```{r}
p <- ggplot(data = sir.output) + 
  geom_line(mapping = aes(x = time, y = S), 
            color = 'blue', size = 2) + 
  geom_line(mapping = aes(x = time, y = R), 
            color = 'black', size = 2) + 
  geom_line(mapping = aes(x = time, y = I), 
            color = 'red', size = 2) + 
  xlim(0,50) + 
  xlab("Time") + 
  ylab("Number of People") + 
  theme(axis.text = element_text(size = 20),
        axis.title = element_text(size = 24)) + 
  theme_bw()

p
```

## Question: What do you notice about the model output? 
Describe what you see about all three state variables. In particular, how does the number of infected people vary over time?

## Question: Varying beta and gamma
Try a few different combinations of beta and gamma, and see how the plot changes.

## Questions: Analyzing the outbreak

How many total people became infected? (The epidemic is over when fewer than one people remain infected)
```{r}
# Find the first time when the number of infected drops below 1
sir.output[sir.output$I <=1,][1,]$R
```

What time was the peak of the epidemic?
```{r}
sir.output[sir.output$I == max(sir.output$I),]$time
```

How many people were infected at the time of the peak?
```{r}
sir.output[sir.output$I == max(sir.output$I),]$I
```