# SIR Model for Covid19

It’s been a hectic few months since my last post. Working and schooling at home during this pandemic chaos is challenging, but it doesn’t stop the hope of renew and transformation. One of the most popular statistcal modeling method for epidemiology and network science-Spreading Phenomena is SIR (susceptible-infected-recovered) model from Tim Churches.

**Susceptible (S)**: Healthy individuals who have not yet contacted the pathogen.

**Infectious (I)**: Contagious individuals who have contacted the pathogen and hence can infect others.

**Recovered/Removed (R)**: Infected individuals become immune or die, i.e., will not be infected again and cannot infect anyone else.

As shown in the picture above, there are two important rates \(\beta\)(transmission rate of the pathogen), \(\gamma\)(recovery rate). You probablily heard of the basic reproduction number \(R_0\) value (R-nought), which can be calculated by \(\frac{\beta}{\gamma}\).

If \(R_0\) > 1, epidemic is in the epidemic state;

If \(R_0\) < 1, epidemic dies out.

**The purpose of the post is to estiamte the basic reproduction number R-nought / \(R_0\) value based on US data and then make predictions for the future.**

**1. Load packages**:

```
#install.packages(c("httr", "jsonlite"))
library(httr)
library(jsonlite)
#install.packages("optimr")
library(optimr)
library(deSolve)
library(tidyverse)
library(lubridate)
```

**2. Retrieve data from John’s Hopkins and use date March 1st till March 15 to fit SIR Model and check if the fit is reasonable**:

```
#by Johns Hopkins CSSE
res = GET("https://pomber.github.io/covid19/timeseries.json")
data = fromJSON(rawToChar(res$content))
#names(data) all contries
US = data$US
US$date=as.Date(US$date)
new<-US %>% dplyr::filter(date>= "2020-03-01"& date<"2020-03-15") #for US data only
head(new) #first six rows
```

```
## date confirmed deaths recovered
## 1 2020-03-01 32 1 7
## 2 2020-03-02 54 6 7
## 3 2020-03-03 74 7 7
## 4 2020-03-04 107 11 7
## 5 2020-03-05 184 12 7
## 6 2020-03-06 237 14 7
```

`tail(new) #last six rows`

```
## date confirmed deaths recovered
## 9 2020-03-09 589 22 7
## 10 2020-03-10 782 28 8
## 11 2020-03-11 1145 33 8
## 12 2020-03-12 1584 43 12
## 13 2020-03-13 2214 51 12
## 14 2020-03-14 2971 58 12
```

**3. build SIR Model and find \(\beta\) and \(\gamma\) by optimization method RSS**:

```
#Initial values:
Infected <- new$confirmed
Day <- 1:(length(new$date))
N<-327200000 #US population
init <- c(S = N - Infected[1], I = Infected[1], R = 0) #Inital value
## SIR Model
SIR <- function(time, state, parameters) {
par <- as.list(c(state, parameters))
with(par, {
dS <- -beta * I * S/N
dI <- beta * I * S/N - gamma * I
dR <- gamma * I
list(c(dS, dI, dR))
})
}
#Optimization by RSS to get beta and gamma
RSS <- function(parameters) {
names(parameters) <- c("beta", "gamma")
out <- ode(y = init, times = Day, func = SIR, parms = parameters)
fit <- out[, 3]
sum((Infected - fit)^2)
}
Opt <- optim(c(0.5, 0.5), RSS, method = "L-BFGS-B", lower = c(0,0), upper = c(1, 1))
Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par
```

```
## beta gamma
## 0.6755947 0.3244053
```

```
#R0
Opt$par[1]/Opt$par[2]
```

`## [1] 2.082563`

**4. Plot the predicted value vs raw data**

```
sir_start_date <- "2020-03-01"
t <- 1:as.integer(ymd("2020-03-15") - ymd(sir_start_date))
# get the fitted values from our SIR model
fitted_cumulative_incidence <- data.frame(ode(y = init, times = t,
func = SIR, parms = Opt_par))
# add a Date column and join the observed incidence data
fitted_cumulative_incidence <- fitted_cumulative_incidence %>%
mutate(date = ymd(sir_start_date) + days(t - 1)) %>%
left_join(new %>% ungroup() %>% select(date, confirmed))
```

`## Joining, by = "date"`

```
# plot the data
fitted_cumulative_incidence %>%
ggplot(aes(x = date)) + geom_line(aes(y = I), colour = "red") +
geom_point(aes(y = confirmed), colour = "orange") +
labs(y = "Cumulative incidence", title = "COVID-19 fitted vs observed cumulative incidence, US 03/01-03/15",
subtitle = "(red=fitted incidence from SIR model, orange=observed incidence)")
```

**5. The fit looks reasonable with R0=1.78, now we use the model to predict the curve for 3 months starting from March**

```
#Prediction
# time in days for predictions
t <- 1:90
# get the fitted values from our SIR model
fitted_cumulative_incidence <- data.frame(ode(y = init, times = t,
func = SIR, parms = Opt_par))
# add a Date column and join the observed incidence data
fitted_cumulative_incidence <- fitted_cumulative_incidence %>%
mutate(date = ymd(sir_start_date) + days(t - 1)) %>%
left_join(new %>% ungroup() %>% select(date, confirmed))
```

`## Joining, by = "date"`

```
# plot the data
fitted_cumulative_incidence %>% ggplot(aes(x = date)) + geom_line(aes(y = I),
colour = "red") + geom_line(aes(y = S), colour = "black") +
geom_line(aes(y = R), colour = "green") + geom_point(aes(y = confirmed),
colour = "orange") + scale_y_continuous(labels = scales::comma) +
labs(y = "Persons", title = "COVID-19 3 months prediction") +
scale_colour_manual(name = "", values = c(red = "red", black = "black",
green = "green", orange = "orange"), labels = c("Susceptible",
"Recovered", "Observed incidence", "Infectious"))
```

**Conclusion: It looks like the peak will be in the end of April and it will die down at the end of May.**