1. Introduction

We are interested in exploring the effect of “social distancing” to “spread of virus”. In particular, we wish to answer the following question: “What if some people do not follow the guideline for social distancing?” Since the society is connected to each other, we think that the effect of disobeying the social distancing might be more severe than we thought.

To check this idea, we make the following assumptions:

  1. The population consists of two groups: One is following the social distancing. The other does not follow the social distancing. We call the former as “normal” group and the latter as “non-normal” group. Let \(W_1\) be the population proportion of the normal group and \(W_2\) be the population proportion of the non-normal group.

  2. A person in the normal group meets \(m_1\) people daily on average. Also, a person in the non-normal group meets \(m_2\) people daily on average, where \(m_2\) is larger than \(m_1\). The actual numbers of the meetings are generated from Poisson distributions. The population average is \(m=W_1 m_1 + W_2 m_2\). We set \(m=4.4\).

  3. If a carrir (= person with virus) meets someone at random, the virus will be tranfered with probability \(p=0.03\).

  4. If someone get the virus then he/she will be available (without self-segregration) for ten day. After 10 days, he/she will be hospotalized.

In the first simulation, we consider three scenarios:

  1. No violation of the social distancing.
  2. Modest violation (20% meets \(m=10\) people daily)
  3. Severe violation (5% meets \(m=31\) people daily)

2. Computation

Notation:

  • \(Q_t\): the number of new cases at the \(t\)-th date.
  • \(C_t\): cumulative number of cases at the \(t\)-th date.
  • \(n_{tj}\): realized number of meetings for group \(j\) at the \(t\)-th date.
  • \(a\): initial cases (we used \(a=10\))
  • \(W_j\): population proportion for group \(j\) (\(j=1\) normal group)
  • \(p_t\): infection rate for each meeting at the \(t\)-th date
  • \(N\): total population size (N=50,000,000)

Formula

We compute the number of cases recursively by the following formula.

\[ Q_t = \left(C_{t-1} - C_{t - 10}\right)\left(\sum_{j = 1}^2 W_j^*n_{tj}p_{(t-1)j}\right) \]

where initial values are

\[ C_{-8} = C_{-7} = \cdots =C_0 = 0 \]

and \(W_j^*\) and \(n_{tj}\) are

\[ W_j^* = \frac{m_jW_j}{\sum_{j= 1}^2{m_jW_j}}, \ \ \ n_{tj} \sim Poi(m_j) \]

Also, the infection probability also changes over time: \[ p_t = p_0 \left(1 - \frac{C_{tj}}{NW_j} \right) =p_0 \left(1 - \frac{C_{t}W_j^*}{NW_j} \right) \]

3. Simulation Study

Simulation Parameters

p_0 = 0.03 # Initial probability
N = 50000000 # Number of population
a = 10 # initial number of cases
B = 1000 # Number of Simulation
set.seed(1)
##### Complex Model ######

m = list(c(4.4),  c(3, 10), c(3, 31))
w = list(c(1),  c(0.8, 0.2), c(0.95, 0.05))
l = 50
res = mapply(resmat, m = m, w = w, MoreArgs = list(l = l, isSimple = FALSE))
yli = min(res[[length(m)]][[1]])
yui = max(res[[length(m)]][[1]])
colm = c(3, 2, 4, 0)
for(i in 1:length(m)){
  boxplot(res[[i]][[1]], outline = FALSE, #main = sprintf("m = %g, w = %g", m[[i]], w[[i]]), 
          col = colm[i], log = "y", ylim = c(yli, yui), ylab = "log-scale cases", xlab = "time")
  Csummary = apply(res[[i]][[1]], 2, function(x){summary(x)[2:5]})
  points(1:l, Csummary["Mean",], type = "l")
  legend("topleft", legend =  sprintf("m = %g, w = %g", m[[i]], w[[i]]))
}

par(mfrow = c(1,1))
for(i in length(m):1){
  boxplot(res[[i]][[1]], outline = FALSE, 
          col = colm[i], add = if(i != length(m)) TRUE else FALSE,
          main = "number of cases w/ different values of m and w", log = "y", ylab = "log-scale cases", xlab = "time")
  Csummary = apply(res[[i]][[1]], 2, function(x){summary(x)[2:5]})
  points(1:l, Csummary["Mean",], type = "l")
}
legend("topleft", legend = paste("m=", m,"w=", w), fill = colm, cex = 1)

Each box plot represents the median, 25%, 75% quantile of the cumulative cases and the solid line represents the mean of simulated results.

# C_t; cumulative number of cases
sapply(sequence(length(m)), function(x){
  apply(res[[x]][[1]], 2, function(x){summary(x)["Mean"]})
})
##            [,1]       [,2]         [,3]
##  [1,]  10.00000   10.00000 1.000000e+01
##  [2,]  11.31070   11.82686 1.384738e+01
##  [3,]  12.78459   14.04300 1.919471e+01
##  [4,]  14.46107   16.63239 2.660583e+01
##  [5,]  16.39074   19.76856 3.694124e+01
##  [6,]  18.58268   23.39759 5.113927e+01
##  [7,]  21.05753   27.75163 7.072836e+01
##  [8,]  23.84220   32.95961 9.808114e+01
##  [9,]  26.97328   39.13852 1.360701e+02
## [10,]  30.51003   46.43635 1.890250e+02
## [11,]  33.30294   53.16629 2.586692e+02
## [12,]  36.14460   60.82196 3.536164e+02
## [13,]  39.20156   69.55682 4.832304e+02
## [14,]  42.56114   79.31733 6.598968e+02
## [15,]  46.01500   90.34822 8.980848e+02
## [16,]  49.68053  102.69649 1.223439e+03
## [17,]  53.49801  116.47445 1.664706e+03
## [18,]  57.39135  132.01026 2.267902e+03
## [19,]  61.44799  149.44922 3.087681e+03
## [20,]  65.47572  168.54763 4.201342e+03
## [21,]  69.68483  190.26540 5.724652e+03
## [22,]  73.98429  214.27419 7.807589e+03
## [23,]  78.57553  241.40642 1.062815e+04
## [24,]  83.36972  271.72105 1.448079e+04
## [25,]  88.27871  305.13703 1.974533e+04
## [26,]  93.61660  342.22979 2.683751e+04
## [27,]  98.96162  383.90946 3.664291e+04
## [28,] 104.34369  430.61570 4.996846e+04
## [29,] 110.14341  482.68726 6.787715e+04
## [30,] 115.93314  541.31009 9.216034e+04
## [31,] 122.05147  606.41340 1.253322e+05
## [32,] 128.48073  679.41956 1.702225e+05
## [33,] 135.05608  761.37113 2.299789e+05
## [34,] 141.81356  853.29013 3.099301e+05
## [35,] 148.85399  954.78907 4.169646e+05
## [36,] 156.02105 1070.27881 5.589242e+05
## [37,] 163.56219 1196.84837 7.449734e+05
## [38,] 171.37307 1339.94645 9.866582e+05
## [39,] 179.35982 1498.81618 1.295026e+06
## [40,] 187.81436 1674.04199 1.682100e+06
## [41,] 196.40267 1868.77226 2.147524e+06
## [42,] 205.60729 2085.41256 2.694916e+06
## [43,] 214.58563 2329.86691 3.313283e+06
## [44,] 224.45517 2607.05196 3.986301e+06
## [45,] 234.29028 2915.09144 4.672108e+06
## [46,] 244.51011 3257.42014 5.330295e+06
## [47,] 255.22960 3638.54939 5.926445e+06
## [48,] 266.16251 4071.51651 6.442207e+06
## [49,] 277.65303 4550.47575 6.861927e+06
## [50,] 289.33603 5082.27833 7.195904e+06

At each row, the expected numbers of total cases are presented for each scenario.

4. Further simulation Study

We also consider other scenarios for the 5% outlier groups:

##### Complex Model ######
m = list(c(3, 30), c(5, 10),  c(3, 20))
w = list(c(0.99, 0.01), c(0.95, 0.05), c(0.95, 0.05))

l = 50
res = mapply(resmat, m = m, w = w, MoreArgs = list(l = l, isSimple = FALSE))
yli = min(res[[length(m)]][[1]])
yui = max(res[[length(m)]][[1]])
colm = c(3, 2, 4, 0)
for(i in 1:length(m)){
  boxplot(res[[i]][[1]], outline = FALSE, #main = sprintf("m = %g, w = %g", m[[i]], w[[i]]), 
          col = colm[i], log = "y", ylim = c(yli, yui), ylab = "log-scale cases", xlab = "time")
  Csummary = apply(res[[i]][[1]], 2, function(x){summary(x)[2:5]})
  points(1:l, Csummary["Mean",], type = "l")
  legend("topleft", legend =  sprintf("m = %g, w = %g", m[[i]], w[[i]]))
}

par(mfrow = c(1,1))
for(i in length(m):1){
  boxplot(res[[i]][[1]], outline = FALSE, 
          col = colm[i], add = if(i != length(m)) TRUE else FALSE,
          main = "number of cases w/ different values of m and w", log = "y", ylab = "log-scale cases", xlab = "time")
  Csummary = apply(res[[i]][[1]], 2, function(x){summary(x)[2:5]})
  points(1:l, Csummary["Mean",], type = "l")
}
legend("topleft", legend = paste("m=", m,"w=", w), fill = colm, cex = 1)

Each box plot represents the median, 25%, 75% quantile of the cumulative cases and the solid line represents the mean of simulated results.

# C_t; cumulative number of cases
sapply(sequence(length(m)), function(x){
  apply(res[[x]][[1]], 2, function(x){summary(x)["Mean"]})
})
##             [,1]       [,2]        [,3]
##  [1,]   10.00000   10.00000    10.00000
##  [2,]   11.64691   11.64154    12.19817
##  [3,]   13.54720   13.52920    14.88076
##  [4,]   15.77245   15.74376    18.21144
##  [5,]   18.31727   18.35275    22.22590
##  [6,]   21.32773   21.38052    27.12712
##  [7,]   24.81730   24.90517    33.12130
##  [8,]   28.87236   29.04948    40.50009
##  [9,]   33.51658   33.70933    49.60788
## [10,]   39.01044   39.29479    60.68660
## [11,]   43.74622   44.17294    72.00406
## [12,]   49.02949   49.53764    85.20096
## [13,]   54.83329   55.41431   100.93925
## [14,]   61.20349   61.97355   119.25801
## [15,]   68.26192   69.12289   140.84144
## [16,]   75.92794   76.95932   166.05353
## [17,]   84.26166   85.51982   195.67483
## [18,]   93.37509   94.92947   230.43701
## [19,]  103.38486  104.84895   270.85360
## [20,]  114.04369  115.69596   317.52200
## [21,]  125.83803  127.49897   372.42601
## [22,]  138.42604  140.32559   436.81195
## [23,]  152.28660  154.54600   511.05325
## [24,]  167.22697  170.28555   597.83866
## [25,]  183.45790  187.12667   700.38844
## [26,]  201.08101  205.17810   820.41803
## [27,]  220.23830  224.74572   958.28555
## [28,]  241.18906  246.29606  1120.39456
## [29,]  264.07701  269.29096  1310.42458
## [30,]  288.55118  294.50668  1532.00151
## [31,]  315.20379  321.90611  1784.94719
## [32,]  344.29742  351.05663  2085.17240
## [33,]  375.80685  382.85678  2434.67460
## [34,]  410.25620  418.14323  2843.22958
## [35,]  447.20437  455.74955  3321.36864
## [36,]  487.12623  497.44043  3874.51319
## [37,]  531.58283  543.08013  4520.30594
## [38,]  579.63483  592.31538  5276.90328
## [39,]  631.32094  645.63831  6153.43431
## [40,]  688.46075  704.76661  7185.76527
## [41,]  748.55005  767.97538  8399.34805
## [42,]  815.21945  836.54510  9807.63314
## [43,]  888.59812  910.13615 11410.83715
## [44,]  967.15267  989.21966 13313.66314
## [45,] 1052.27655 1074.09353 15543.97847
## [46,] 1143.99207 1170.60967 18126.00004
## [47,] 1244.89918 1272.06464 21169.67365
## [48,] 1354.01844 1385.62421 24711.38685
## [49,] 1470.54379 1505.74482 28875.41396
## [50,] 1598.56997 1640.82700 33624.23560

5. Discussion

Findings (from first simulation)

  1. Under scenario 1 (everybody follows the social distancing), the total cases increase only 29 times after 50 days.
  2. Under scenario 2 (80% follows social distancing, but 20% make modest violation by meeting 10 people every day), the total cases increase 508 times after 50 days.
  3. Under scenario 3 (95% follows social distancing, but 5% make severe violation by meeting 31 people every day), the total cases increase more than 719,590 times after 50 days.
  4. On average, the above three scenarios have the same value of average meetings (\(m=4.4\)), but the existence of the outlying groups makes a big difference in the total number of cases.

Findings (from the second simulation)

  1. Under Scenario 1 (95% follows \(m=5\) and 5% follows \(m=10\)), the total number of cases increases to 159 times.
  2. Under scenorio 2 (99% follows \(m=3\) and 1% follows \(m=30\)), the total number of cases is about the same as the first scenario.
  3. Under Scenario 3 (95% follows \(m=3\) and 5% follows \(m=20\)), the total number of cases increases to 3362 times.

For comments or questions, please contact "jkim@iastate.edu"

Social distancing is critical!