#Packages
library(dplyr); library(magrittr); 
library(ggplot2); library(ggsci)
library(reshape2)
theme_set(theme_minimal())

This document serves as a technical companion to the piece titled Flattning the Curve at OR/MS Today. As a starter, we repeat the Kermack-McKendric S-I-R model,

\[ dS/dt = -\alpha SI \\ dI/dt = \alpha SI -\rho I \\ dR/dt = \rho I .\]

This is a classic ‘cohort model’; written in matrix form, we can see that the system is closed and the total differential (sum over columns when properly formatted) is zero.

Of interest are the parameters \(\alpha\), which is a measure of both a virus’ infectivity and speed of spread, and \(\rho\) which is a measure (Markovian) of how long the infectious period lasts.

Here’s a function that returns a time-trace of the SIR Model from given Initial Conditions. h is a parameter that controls the internal time-step; i.e. \(dx/dt = lim \space h \rightarrow 0 \space (x(t+h)-x(t))/h\) We are using Euler’s method to evaluate this system of ordinary differential equations. Smaller values of h lead to more accuracy (and also more computing time!)

KMC = function(So = .99, Io = .01, Ro = 0, h = .1, 
               alpha = .1, rho = 1/14, tmax = 365){
  Time = seq(0, tmax, by = h)
  #Normalize
  Tot = So + Io + Ro
  S = So/Tot; I = Io/Tot ;R = Ro/Tot
  St = vector(); It = vector(); Rt = vector()
  St[1] = S; It[1] = I; Rt[1] = R
  #Loop
  for(i in 2:length(Time)){
    St[i] = max(0, St[i-1] - 
                  (St[i-1]* It[i-1]*alpha *h))
    It[i] = max(0, It[i-1] + 
                  (St[i-1]* It[i-1]*alpha *h) - (rho*h*It[i-1]))
    Rt[i] = Rt[i-1] + 
      (rho*h*It[i-1])
  }
  res = data.frame(t = Time, S = St, I = It, R = Rt)
  return(res) 
}

Below is a plot of the process, accepting all default values from above:

X = KMC()

X %>% melt(id.vars = "t") %>% ggplot(aes(x = t, y = value, color = variable)) + 
  geom_line() + scale_color_manual(values = c("blue", "red", "green")) + 
  xlab("Time (days)") + scale_y_continuous(labels = scales::percent, name = "Proportion of Population" )  + labs(color = "Population")

Parameter Estimation

Our approach is as follows: We wish to vary the virulance of the disease, \(\alpha\) against a fixed, realistic choice for Removals, \(R\). Thus, we wish to benchmark the sojourn time of any given individual to the disease in the presumed incubation period, which is currently estimated at 14 days (Carol check this). While there are several ways to fit data to an epidemic, we will propose the following straightforward application of differential equations.

Fitting \(\rho\)

Let’s consider the equation for I(t):

\[ dI/dt = \alpha SI - \rho I \]

Imagine, if you will, that you have just been infected, and you have immediately been placed alone in quarantine until you are no longer infective. Therefore, the equation above collapses to:

\[ dI/dt = -\rho I \\ s.t. \\ I(0) = 1 \]

Which has the solution

\[I(t) = 1- e^{-\rho t}, \] which the probabilisits will immediately recognize as the CDF of the exponential distribution (this is not a coincidence!) and given estimates of a 2-week infective period 1, we conclude that \(\rho \approx .07\)

Estimating \(\alpha\)

We use data from the CDC (table at bottom of article),

CD = read.csv("CDCData.csv")

CD %>% ggplot(aes(x = Julian, y = Cumulative)) + geom_line() + xlab("Julian Date 2020") + ylab("Cumulative cases") 

CDD = CD %>% filter(Julian > 45) #Truncates early effects

CDD$Log = log(CDD$Cumulative)
CDD$JD = CDD$Julian - 45

lm1 = lm(Log~JD, data = CDD)

lm1 %>% broom::tidy() %>% knitr::kable()
term estimate std.error statistic p.value
(Intercept) 3.1072382 0.0759074 40.93459 0
JD 0.1036235 0.0055415 18.69961 0

We will use \(\alpha = .1\) as a starting point for analysis. This is consistent with a number distilled separately by a group of researchers at LANL 2

Note that the LANL paper uses \(\alpha = .1, R_0 \approx 2.3\) based on infections in Wuhan, China after control measures, with an implied \(\rho = .04\) Our estimate is roughly

Cases on \(\alpha\)

The point of this is not the specific values of \(\alpha\) or \(\rho\), but rather how incremental reductions in \(\alpha\) can change the

rho = 1/14;alpha_0 = .1; tmax = 2*365

effects = c(0, .01, .05, .07,.1, .15, .25, .35)
names = c("t", paste0(100*effects,"Pct"))
names = c("t", paste0(100*effects,"%"))
n = length(effects)

X = KMC(tmax = tmax)

output = data.frame(t = X$t) #Infection Curve
out2 = vector() #Total Infection amount

for(i in 1:n){
  X = KMC(So = .9999, Io = 1-.9999, rho = rho, alpha = alpha_0*(1-effects[i]), tmax = tmax)
  output %<>% cbind(X$I)
  out2[i] = max(X$R)
}

names(output) = names

output %>% melt(id.vars = "t") %>% 
  ggplot(aes(x = t, y = value, color = variable)) + geom_line() + 
  xlab("time (days)") + scale_y_continuous(labels = scales::percent, name = "Infected") + labs(color = "Countermeasure Effectiveness")

Ultimate size of the epidemic:

data.frame(effect = effects, Size = out2) %>% ggplot(aes(x = effect, y= Size)) + geom_line() + scale_y_continuous(name = "Ultimate size of Epidemic", labels = scales::percent) + xlab("Effectiveness of Countermeasures") + scale_x_continuous(labels = scales::percent)

Data from CDC:

library(knitr)
CD %>% kable()
Julian New.Cases Cumulative
11 0 0
12 0 0
13 2 2
14 0 2
15 1 3
16 0 3
17 0 3
18 0 3
19 0 3
20 1 4
21 1 5
22 0 5
23 1 6
24 2 8
25 0 8
26 0 8
27 2 10
28 2 12
29 0 12
30 0 12
30 0 12
31 0 12
32 0 12
33 0 12
34 0 12
35 0 12
36 0 12
37 0 12
38 0 12
39 1 13
40 1 14
41 0 14
42 2 16
43 1 17
44 0 17
45 0 17
46 3 20
47 4 24
48 1 25
49 3 28
50 6 34
51 9 43
52 5 48
53 11 59
54 5 64
55 14 78
56 6 84
57 13 97
58 9 106
60 19 125
61 10 135
62 9 144
63 7 151
64 12 163
65 2 165
66 7 172
67 0 172
68 0 172

  1. “we have little confidence in these estimates”↩

  2. “https://www.medrxiv.org/content/10.1101/2020.02.07.20021154v1.full.pdf”↩