Simple SIR data simulation
Breck Baldwin
May 18, 2021
Model 1
About the simplest SIR model I can think of. Surprisingly error prone process and very sensitive to infection rate. See below. I find it very useful to naively code up a simulation for some problem I am trying to model. I am not an epidemiologist and have no medical training.
I have looked at some infection models so this is not totally naive, see https://mc-stan.org/users/documentation/case-studies/boarding_school_case_study.html which I found to be a good explanation. I have not attempted to reproduce that model here although there is a model reproduction check list I did at https://codatmo.github.io/Simple_SIR/.
I donāt use standard terms like R0, Rt etc because there are subtleties that I am not sure will get right and wanting to keep things naive.
There are likely gross errors which please point out in the issue tracker or send me email.
rm(list = ls()) #clean out the cruft
set.seed(43614)
runModel1 = function(runName = runName,
nPop = nPop,
nDays = nDays,
print = PRINT,
probInfectPerContact = probInfectPerContact) {
dayState = c('i', rep('s', nPop - 1))
infectionDay = c(1, rep(NA, nPop - 1))
meanDaysInfectious = 3
meanContactsPerDay = 10
#dataM = matrix(nrow = nDays, ncol = 3)
colNames = c('runName', 'day', 's', 'i', 'r', 'probInfectPerContact')
df = data.frame(matrix(nrow = 0, ncol=length(colNames)))
colnames(df) = colNames
nextDayState = dayState
for (d in 1:nDays) {
df[d,] = rep(NA,length(colNames)) #setup data, will cause errors if I miss something
df[d,]$runName = runName
df[d,]$s= length(subset(dayState, dayState == 's'))
df[d,]$i = length(subset(dayState, dayState == 'i'))
df[d,]$r = length(subset(dayState, dayState == 'r'))
df[d,]$day = d
df[d,]$probInfectPerContact = probInfectPerContact
#dataM[d,] = c(susceptible, infectious, resolved)
if (print) {
cat(
sprintf(
"Day=%d, susceptible=%d, infected=%d, resolved=%d\n",
df[d,]$day, df[d,]$s, df[d,]$i, df[d,]$r))
}
for (per in 1:nPop) {
#end infectious period
if (dayState[per] == 'i') {
if (d - infectionDay[per] > meanDaysInfectious) {
nextDayState[per] = 'r'
}
}
}
for (per in 1:nPop) {
if (dayState[per] == 'i') {
for (otherPer in sample(1:nPop, meanContactsPerDay)) {
if (dayState[otherPer] == 's' &&
rbinom(1, 1, probInfectPerContact) == 1) {
nextDayState[otherPer] = 'i'
infectionDay[otherPer] = d
}
}
}
}
dayState = nextDayState
}
return(df)
}
dataM = runModel1(runName = "1/30", nPop = 1e+04, nDays = 100, print = TRUE, probInfectPerContact = 1/30)Day=1, susceptible=9999, infected=1, resolved=0
Day=2, susceptible=9998, infected=2, resolved=0
Day=3, susceptible=9997, infected=3, resolved=0
Day=4, susceptible=9996, infected=4, resolved=0
Day=5, susceptible=9996, infected=4, resolved=0
Day=6, susceptible=9995, infected=3, resolved=2
Day=7, susceptible=9993, infected=4, resolved=3
Day=8, susceptible=9993, infected=3, resolved=4
Day=9, susceptible=9991, infected=5, resolved=4
Day=10, susceptible=9989, infected=6, resolved=5
Day=11, susceptible=9989, infected=4, resolved=7
Day=12, susceptible=9989, infected=4, resolved=7
Day=13, susceptible=9988, infected=3, resolved=9
Day=14, susceptible=9987, infected=2, resolved=11
Day=15, susceptible=9986, infected=3, resolved=11
Day=16, susceptible=9984, infected=5, resolved=11
Day=17, susceptible=9983, infected=5, resolved=12
Day=18, susceptible=9979, infected=8, resolved=13
Day=19, susceptible=9974, infected=12, resolved=14
Day=20, susceptible=9971, infected=13, resolved=16
Day=21, susceptible=9965, infected=18, resolved=17
Day=22, susceptible=9959, infected=20, resolved=21
Day=23, susceptible=9949, infected=25, resolved=26
Day=24, susceptible=9946, infected=25, resolved=29
Day=25, susceptible=9942, infected=23, resolved=35
Day=26, susceptible=9932, infected=27, resolved=41
Day=27, susceptible=9923, infected=26, resolved=51
Day=28, susceptible=9917, infected=29, resolved=54
Day=29, susceptible=9906, infected=36, resolved=58
Day=30, susceptible=9893, infected=39, resolved=68
Day=31, susceptible=9878, infected=45, resolved=77
Day=32, susceptible=9863, infected=54, resolved=83
Day=33, susceptible=9841, infected=65, resolved=94
Day=34, susceptible=9818, infected=75, resolved=107
Day=35, susceptible=9793, infected=85, resolved=122
Day=36, susceptible=9762, infected=101, resolved=137
Day=37, susceptible=9730, infected=111, resolved=159
Day=38, susceptible=9694, infected=124, resolved=182
Day=39, susceptible=9648, infected=145, resolved=207
Day=40, susceptible=9604, infected=158, resolved=238
Day=41, susceptible=9563, infected=167, resolved=270
Day=42, susceptible=9519, infected=175, resolved=306
Day=43, susceptible=9472, infected=176, resolved=352
Day=44, susceptible=9424, infected=180, resolved=396
Day=45, susceptible=9365, infected=198, resolved=437
Day=46, susceptible=9305, infected=214, resolved=481
Day=47, susceptible=9251, infected=221, resolved=528
Day=48, susceptible=9183, infected=241, resolved=576
Day=49, susceptible=9118, infected=247, resolved=635
Day=50, susceptible=9044, infected=261, resolved=695
Day=51, susceptible=8967, infected=284, resolved=749
Day=52, susceptible=8879, infected=304, resolved=817
Day=53, susceptible=8806, infected=312, resolved=882
Day=54, susceptible=8717, infected=327, resolved=956
Day=55, susceptible=8628, infected=339, resolved=1033
Day=56, susceptible=8541, infected=338, resolved=1121
Day=57, susceptible=8439, infected=367, resolved=1194
Day=58, susceptible=8333, infected=384, resolved=1283
Day=59, susceptible=8232, infected=396, resolved=1372
Day=60, susceptible=8122, infected=419, resolved=1459
Day=61, susceptible=7994, infected=445, resolved=1561
Day=62, susceptible=7896, infected=437, resolved=1667
Day=63, susceptible=7779, infected=453, resolved=1768
Day=64, susceptible=7674, infected=448, resolved=1878
Day=65, susceptible=7576, infected=418, resolved=2006
Day=66, susceptible=7462, infected=434, resolved=2104
Day=67, susceptible=7360, infected=419, resolved=2221
Day=68, susceptible=7247, infected=427, resolved=2326
Day=69, susceptible=7135, infected=441, resolved=2424
Day=70, susceptible=7044, infected=418, resolved=2538
Day=71, susceptible=6955, infected=405, resolved=2640
Day=72, susceptible=6862, infected=385, resolved=2753
Day=73, susceptible=6784, infected=351, resolved=2865
Day=74, susceptible=6720, infected=324, resolved=2956
Day=75, susceptible=6649, infected=306, resolved=3045
Day=76, susceptible=6595, infected=267, resolved=3138
Day=77, susceptible=6528, infected=256, resolved=3216
Day=78, susceptible=6472, infected=248, resolved=3280
Day=79, susceptible=6421, infected=228, resolved=3351
Day=80, susceptible=6370, infected=225, resolved=3405
Day=81, susceptible=6319, infected=209, resolved=3472
Day=82, susceptible=6266, infected=206, resolved=3528
Day=83, susceptible=6222, infected=199, resolved=3579
Day=84, susceptible=6176, infected=194, resolved=3630
Day=85, susceptible=6138, infected=181, resolved=3681
Day=86, susceptible=6100, infected=166, resolved=3734
Day=87, susceptible=6071, infected=151, resolved=3778
Day=88, susceptible=6040, infected=136, resolved=3824
Day=89, susceptible=6008, infected=130, resolved=3862
Day=90, susceptible=5971, infected=129, resolved=3900
Day=91, susceptible=5941, infected=130, resolved=3929
Day=92, susceptible=5923, infected=117, resolved=3960
Day=93, susceptible=5891, infected=117, resolved=3992
Day=94, susceptible=5861, infected=110, resolved=4029
Day=95, susceptible=5834, infected=107, resolved=4059
Day=96, susceptible=5818, infected=105, resolved=4077
Day=97, susceptible=5797, infected=94, resolved=4109
Day=98, susceptible=5775, infected=86, resolved=4139
Day=99, susceptible=5756, infected=78, resolved=4166
Day=100, susceptible=5734, infected=84, resolved=4182Running this showed lots of sensitivity around 1/31 to 1/33 infection rate.
Thought I should graph it, see below, this is for the āiā state or infectious.
library(tidyverse)
library(ggplot2)
df = data.frame()
for (i in 20:34) {
df2 = runModel1(runName = sprintf("1/%d=%.3f", i, 1/i), nPop = 1e+04,
nDays = 100, print = FALSE, probInfectPerContact = 1/i)
df = rbind(df,df2)
}
dfLong = gather(df, key = sir, value = count, c(s, i, r))
plot = ggplot(data = dfLong[dfLong$sir == 'i',],
aes(x = day, y = count, group = runName, color = runName)) +
geom_line()
print(plot)
Above are varying rates of probability of infection per contact for a population of 10,000 for 100 days for the āiā or infected state. Whether patient-zero manages to infect someone is a big issue, real world probably like this too.
nDays = 100
nPop = 1e+04
dataM = runModel1(runName = "1/30", nPop = nPop, nDays = nDays, print = TRUE,
probInfectPerContact = 1/30)Day=1, susceptible=9999, infected=1, resolved=0
Day=2, susceptible=9999, infected=1, resolved=0
Day=3, susceptible=9999, infected=1, resolved=0
Day=4, susceptible=9998, infected=2, resolved=0
Day=5, susceptible=9998, infected=2, resolved=0
Day=6, susceptible=9995, infected=4, resolved=1
Day=7, susceptible=9992, infected=7, resolved=1
Day=8, susceptible=9990, infected=8, resolved=2
Day=9, susceptible=9990, infected=8, resolved=2
Day=10, susceptible=9988, infected=7, resolved=5
Day=11, susceptible=9983, infected=9, resolved=8
Day=12, susceptible=9978, infected=12, resolved=10
Day=13, susceptible=9972, infected=18, resolved=10
Day=14, susceptible=9965, infected=23, resolved=12
Day=15, susceptible=9952, infected=31, resolved=17
Day=16, susceptible=9947, infected=31, resolved=22
Day=17, susceptible=9934, infected=38, resolved=28
Day=18, susceptible=9926, infected=39, resolved=35
Day=19, susceptible=9913, infected=39, resolved=48
Day=20, susceptible=9903, infected=44, resolved=53
Day=21, susceptible=9888, infected=46, resolved=66
Day=22, susceptible=9872, infected=54, resolved=74
Day=23, susceptible=9852, infected=61, resolved=87
Day=24, susceptible=9836, infected=67, resolved=97
Day=25, susceptible=9816, infected=72, resolved=112
Day=26, susceptible=9785, infected=87, resolved=128
Day=27, susceptible=9749, infected=103, resolved=148
Day=28, susceptible=9717, infected=119, resolved=164
Day=29, susceptible=9677, infected=139, resolved=184
Day=30, susceptible=9628, infected=157, resolved=215
Day=31, susceptible=9574, infected=175, resolved=251
Day=32, susceptible=9526, infected=191, resolved=283
Day=33, susceptible=9476, infected=201, resolved=323
Day=34, susceptible=9419, infected=209, resolved=372
Day=35, susceptible=9349, infected=225, resolved=426
Day=36, susceptible=9291, infected=235, resolved=474
Day=37, susceptible=9230, infected=246, resolved=524
Day=38, susceptible=9154, infected=265, resolved=581
Day=39, susceptible=9071, infected=278, resolved=651
Day=40, susceptible=8992, infected=299, resolved=709
Day=41, susceptible=8907, infected=323, resolved=770
Day=42, susceptible=8794, infected=360, resolved=846
Day=43, susceptible=8690, infected=381, resolved=929
Day=44, susceptible=8594, infected=398, resolved=1008
Day=45, susceptible=8461, infected=446, resolved=1093
Day=46, susceptible=8326, infected=468, resolved=1206
Day=47, susceptible=8187, infected=503, resolved=1310
Day=48, susceptible=8045, infected=549, resolved=1406
Day=49, susceptible=7897, infected=564, resolved=1539
Day=50, susceptible=7739, infected=587, resolved=1674
Day=51, susceptible=7602, infected=585, resolved=1813
Day=52, susceptible=7464, infected=581, resolved=1955
Day=53, susceptible=7331, infected=566, resolved=2103
Day=54, susceptible=7179, infected=560, resolved=2261
Day=55, susceptible=7061, infected=541, resolved=2398
Day=56, susceptible=6954, infected=510, resolved=2536
Day=57, susceptible=6833, infected=498, resolved=2669
Day=58, susceptible=6728, infected=451, resolved=2821
Day=59, susceptible=6624, infected=437, resolved=2939
Day=60, susceptible=6522, infected=432, resolved=3046
Day=61, susceptible=6431, infected=402, resolved=3167
Day=62, susceptible=6340, infected=388, resolved=3272
Day=63, susceptible=6263, infected=361, resolved=3376
Day=64, susceptible=6192, infected=330, resolved=3478
Day=65, susceptible=6088, infected=343, resolved=3569
Day=66, susceptible=6030, infected=310, resolved=3660
Day=67, susceptible=5970, infected=293, resolved=3737
Day=68, susceptible=5903, infected=289, resolved=3808
Day=69, susceptible=5831, infected=257, resolved=3912
Day=70, susceptible=5781, infected=249, resolved=3970
Day=71, susceptible=5737, infected=233, resolved=4030
Day=72, susceptible=5688, infected=215, resolved=4097
Day=73, susceptible=5654, infected=177, resolved=4169
Day=74, susceptible=5624, infected=157, resolved=4219
Day=75, susceptible=5599, infected=138, resolved=4263
Day=76, susceptible=5575, infected=113, resolved=4312
Day=77, susceptible=5551, infected=103, resolved=4346
Day=78, susceptible=5532, infected=92, resolved=4376
Day=79, susceptible=5515, infected=84, resolved=4401
Day=80, susceptible=5500, infected=75, resolved=4425
Day=81, susceptible=5482, infected=69, resolved=4449
Day=82, susceptible=5473, infected=59, resolved=4468
Day=83, susceptible=5463, infected=52, resolved=4485
Day=84, susceptible=5451, infected=49, resolved=4500
Day=85, susceptible=5445, infected=37, resolved=4518
Day=86, susceptible=5437, infected=36, resolved=4527
Day=87, susceptible=5430, infected=33, resolved=4537
Day=88, susceptible=5424, infected=27, resolved=4549
Day=89, susceptible=5422, infected=23, resolved=4555
Day=90, susceptible=5419, infected=18, resolved=4563
Day=91, susceptible=5411, infected=19, resolved=4570
Day=92, susceptible=5409, infected=15, resolved=4576
Day=93, susceptible=5406, infected=16, resolved=4578
Day=94, susceptible=5405, infected=14, resolved=4581
Day=95, susceptible=5401, infected=10, resolved=4589
Day=96, susceptible=5400, infected=9, resolved=4591
Day=97, susceptible=5399, infected=7, resolved=4594
Day=98, susceptible=5397, infected=8, resolved=4595
Day=99, susceptible=5396, infected=5, resolved=4599
Day=100, susceptible=5396, infected=4, resolved=4600library("cmdstanr")
# Run from command line: Rscript run.R
# If running from RStudio remember to set the working directory
# >Session>Set Working Directory>To Source File Location
model <- cmdstan_model("stan/sir_negbin.stan")
stan_data <- list(n_days = nrow(dataM), y0 = c(nPop -1, 1, 0), t0 = 0,
ts = 1:nDays, N = nPop, cases = dataM$i, compute_likelihood = 1)
fit <- model$sample(data=stan_data, output_dir="output", num_cores = 4,
num_chains = 4)Warning in model$sample(data = stan_data, output_dir = "output", num_cores =
4, : 'num_cores' is deprecated. Please use 'parallel_chains' instead.Warning in model$sample(data = stan_data, output_dir = "output", num_cores =
4, : 'num_chains' is deprecated. Please use 'chains' instead.Running MCMC with 4 parallel chains...
Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup) Warning: Chain 1 finished unexpectedly!Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2 finished in 16.2 seconds.
Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3 finished in 20.4 seconds.
Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4 finished in 24.6 seconds.Warning: 1 chain(s) finished unexpectedly!The remaining chains had a mean execution time of 30.4 seconds.Warning: The returned fit object will only read in results of successful
chains. Please use read_cmdstan_csv() to read the results of the failed chains
separately.Use the $output(chain_id) method for more output of the failed
chains.print(fit$summary(variables = c('pred_cases')))# A tibble: 100 x 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 pred_cases[1] 1.20 1 1.13 1.48 0 3 1.00 2909. 2859.
2 pred_cases[2] 1.38 1 1.24 1.48 0 4 1.00 2902. 2821.
3 pred_cases[3] 1.63 1 1.37 1.48 0 4 1.00 2904. 2834.
4 pred_cases[4] 1.92 2 1.49 1.48 0 5 1.00 2501. 2775.
5 pred_cases[5] 2.30 2 1.63 1.48 0 5 1.00 2910. 2814.
6 pred_cases[6] 2.67 2 1.77 1.48 0 6 1.00 3104. 3025.
7 pred_cases[7] 3.15 3 1.97 1.48 0 7 1.00 2959. 2507.
8 pred_cases[8] 3.75 3 2.22 1.48 1 8 1.00 2983. 2753.
9 pred_cases[9] 4.36 4 2.42 2.97 1 9 1.00 3055. 2929.
10 pred_cases[10] 5.13 5 2.66 2.97 1 10 1.00 3049. 2995.
# ā¦ with 90 more rowsAdding tweets
Below is a very simple data generating process that produces tweets in addition to modeling a SIR compartment model.
rm(list = ls())
set.seed(43614)
runSim2 = function(runName = runName,
nPop = nPop,
nDays = nDays,
print = PRINT,
probInfectPerContact = probInfectPerContact,
probTweet = probTweet,
probDeath = probDeath) {
dayState = c('i', rep('s', nPop - 1))
infectionDay = c(1, rep(NA, nPop - 1))
tweets = rep(0, nPop)
deaths = rep(0,nPop)
meanDaysInfectious = 3
meanContactsPerDay = 10
#dataM = matrix(nrow = nDays, ncol = 3)
colNames = c('runName', 'day', 's', 'i', 'r', 'probInfectPerContact', 'tweets','deaths')
df = data.frame(matrix(nrow = 0, ncol=length(colNames)))
colnames(df) = colNames
nextDayState = dayState
for (d in 1:nDays) {
df[d,] = rep(NA,length(colNames)) #setup data, will cause errors if I miss something
df[d,]$runName = runName
df[d,]$s= length(subset(dayState, dayState == 's'))
df[d,]$i = length(subset(dayState, dayState == 'i'))
df[d,]$r = length(subset(dayState, dayState == 'r'))
df[d,]$day = d
df[d,]$probInfectPerContact = probInfectPerContact
df[d,]$tweets = sum(tweets)
df[d,]$deaths = sum(deaths)
#dataM[d,] = c(susceptible, infectious, resolved)
if (print) {
cat(
sprintf(
"Day=%d, susceptible=%d, infected=%d, resolved=%d, tweets=%d, deaths=%d\n",
df[d,]$day, df[d,]$s, df[d,]$i, df[d,]$r, df[d,]$tweets, df[d,]$deaths))
}
tweets = rep(0, nPop) #start fresh every day, certainly wrong.
for (per in 1:nPop) {
#end infectious period
if (dayState[per] == 'i') {
if (d - infectionDay[per] > meanDaysInfectious) {
nextDayState[per] = 'r'
if (rbinom(1,1,probDeath) == 1) { #do they die?
deaths[per] = 1;
}
}
}
}
for (per in 1:nPop) {
if (dayState[per] == 'i') {
tweets[per] = rbinom(1, 1, probTweet)
for (otherPer in sample(1:nPop, meanContactsPerDay)) {
if (dayState[otherPer] == 's' &&
rbinom(1, 1, probInfectPerContact) == 1) {
nextDayState[otherPer] = 'i'
infectionDay[otherPer] = d
}
}
}
}
dayState = nextDayState
}
return(df)
}
nWeeks = 10
nDays = nWeeks * 7
nPop = 1e+04
deathsTweetsDf = runSim2(runName = "1/30", nPop = nPop, nDays = nDays, print = TRUE,
probInfectPerContact = 1/30,
probTweet = 1/10,
probDeath = 3/100)Day=1, susceptible=9999, infected=1, resolved=0, tweets=0, deaths=0
Day=2, susceptible=9998, infected=2, resolved=0, tweets=0, deaths=0
Day=3, susceptible=9997, infected=3, resolved=0, tweets=0, deaths=0
Day=4, susceptible=9996, infected=4, resolved=0, tweets=1, deaths=0
Day=5, susceptible=9996, infected=4, resolved=0, tweets=3, deaths=0
Day=6, susceptible=9994, infected=4, resolved=2, tweets=0, deaths=0
Day=7, susceptible=9993, infected=4, resolved=3, tweets=0, deaths=0
Day=8, susceptible=9992, infected=4, resolved=4, tweets=2, deaths=0
Day=9, susceptible=9988, infected=8, resolved=4, tweets=1, deaths=0
Day=10, susceptible=9986, infected=8, resolved=6, tweets=0, deaths=0
Day=11, susceptible=9985, infected=8, resolved=7, tweets=0, deaths=0
Day=12, susceptible=9982, infected=10, resolved=8, tweets=1, deaths=0
Day=13, susceptible=9977, infected=11, resolved=12, tweets=1, deaths=0
Day=14, susceptible=9972, infected=14, resolved=14, tweets=2, deaths=0
Day=15, susceptible=9966, infected=19, resolved=15, tweets=2, deaths=0
Day=16, susceptible=9957, infected=25, resolved=18, tweets=1, deaths=0
Day=17, susceptible=9943, infected=34, resolved=23, tweets=4, deaths=0
Day=18, susceptible=9932, infected=40, resolved=28, tweets=4, deaths=1
Day=19, susceptible=9917, infected=49, resolved=34, tweets=1, deaths=2
Day=20, susceptible=9899, infected=58, resolved=43, tweets=3, deaths=3
Day=21, susceptible=9879, infected=64, resolved=57, tweets=4, deaths=4
Day=22, susceptible=9860, infected=72, resolved=68, tweets=6, deaths=5
Day=23, susceptible=9841, infected=76, resolved=83, tweets=5, deaths=6
Day=24, susceptible=9818, infected=81, resolved=101, tweets=7, deaths=6
Day=25, susceptible=9790, infected=89, resolved=121, tweets=11, deaths=6
Day=26, susceptible=9761, infected=99, resolved=140, tweets=8, deaths=8
Day=27, susceptible=9732, infected=109, resolved=159, tweets=9, deaths=9
Day=28, susceptible=9689, infected=129, resolved=182, tweets=6, deaths=9
Day=29, susceptible=9650, infected=140, resolved=210, tweets=14, deaths=9
Day=30, susceptible=9607, infected=154, resolved=239, tweets=13, deaths=10
Day=31, susceptible=9549, infected=183, resolved=268, tweets=15, deaths=10
Day=32, susceptible=9484, infected=205, resolved=311, tweets=28, deaths=12
Day=33, susceptible=9417, infected=233, resolved=350, tweets=23, deaths=13
Day=34, susceptible=9341, infected=266, resolved=393, tweets=20, deaths=15
Day=35, susceptible=9256, infected=293, resolved=451, tweets=30, deaths=15
Day=36, susceptible=9176, infected=308, resolved=516, tweets=41, deaths=17
Day=37, susceptible=9083, infected=334, resolved=583, tweets=27, deaths=19
Day=38, susceptible=9005, infected=336, resolved=659, tweets=37, deaths=21
Day=39, susceptible=8910, infected=346, resolved=744, tweets=28, deaths=22
Day=40, susceptible=8822, infected=354, resolved=824, tweets=39, deaths=23
Day=41, susceptible=8728, infected=355, resolved=917, tweets=35, deaths=27
Day=42, susceptible=8632, infected=373, resolved=995, tweets=46, deaths=29
Day=43, susceptible=8532, infected=378, resolved=1090, tweets=40, deaths=31
Day=44, susceptible=8419, infected=403, resolved=1178, tweets=31, deaths=33
Day=45, susceptible=8299, infected=429, resolved=1272, tweets=42, deaths=36
Day=46, susceptible=8188, infected=444, resolved=1368, tweets=39, deaths=37
Day=47, susceptible=8055, infected=477, resolved=1468, tweets=50, deaths=40
Day=48, susceptible=7942, infected=477, resolved=1581, tweets=38, deaths=44
Day=49, susceptible=7827, infected=472, resolved=1701, tweets=62, deaths=48
Day=50, susceptible=7719, infected=469, resolved=1812, tweets=46, deaths=50
Day=51, susceptible=7592, infected=463, resolved=1945, tweets=45, deaths=53
Day=52, susceptible=7473, infected=469, resolved=2058, tweets=50, deaths=57
Day=53, susceptible=7352, infected=475, resolved=2173, tweets=50, deaths=61
Day=54, susceptible=7226, infected=493, resolved=2281, tweets=49, deaths=63
Day=55, susceptible=7105, infected=487, resolved=2408, tweets=50, deaths=65
Day=56, susceptible=6975, infected=498, resolved=2527, tweets=39, deaths=71
Day=57, susceptible=6883, infected=469, resolved=2648, tweets=47, deaths=73
Day=58, susceptible=6789, infected=437, resolved=2774, tweets=46, deaths=75
Day=59, susceptible=6692, infected=413, resolved=2895, tweets=43, deaths=77
Day=60, susceptible=6603, infected=372, resolved=3025, tweets=42, deaths=82
Day=61, susceptible=6526, infected=357, resolved=3117, tweets=36, deaths=84
Day=62, susceptible=6431, infected=358, resolved=3211, tweets=30, deaths=86
Day=63, susceptible=6346, infected=346, resolved=3308, tweets=32, deaths=89
Day=64, susceptible=6277, infected=326, resolved=3397, tweets=35, deaths=92
Day=65, susceptible=6205, infected=321, resolved=3474, tweets=25, deaths=95
Day=66, susceptible=6153, infected=278, resolved=3569, tweets=38, deaths=97
Day=67, susceptible=6095, infected=251, resolved=3654, tweets=34, deaths=100
Day=68, susceptible=6038, infected=239, resolved=3723, tweets=24, deaths=102
Day=69, susceptible=5981, infected=224, resolved=3795, tweets=29, deaths=104
Day=70, susceptible=5936, infected=217, resolved=3847, tweets=25, deaths=108deaths = colSums(matrix(deathsTweetsDf$deaths, nrow=7))
twitter_symptons <- deathsTweetsDf$tweetsModel fitting, NOT WORKING
Below is the setup for the Liverpool/Brazil model, I still have some issues with beta_right being out of range.
T <- nDays
maxTime <- T
initial_time <- 0
n_disease_states <- 8
n_beta_pieces <- 20
twitter_symptons_start <- 30
n_rho_twitter_symptons_pieces <- 10
population <- nPop
beta_pieceLength <- ceiling(maxTime / n_beta_pieces)
beta_left_t <- seq(0, maxTime - 1, by = beta_pieceLength)
beta_right_t <- c(beta_left_t[2:n_beta_pieces], maxTime + 1)
real_data_length = length(beta_left_t) + length(beta_right_t) + 1
real_data = c(beta_left_t, beta_right_t, population)
integer_data = c(
maxTime,
length(beta_left_t),
length(beta_right_t),
length(beta_left_t),
n_disease_states
)
integer_data_length <- length(integer_data)
rho_calls_pieceLength <-
ceiling((maxTime - twitter_symptons_start) / n_rho_twitter_symptons_pieces)
rho_twitter_symptons_left_t <-
seq(twitter_symptons_start, maxTime - 1, by = rho_calls_pieceLength)
rho_twitter_symptons_right_t <-
c(rho_twitter_symptons_left_t[2:n_rho_twitter_symptons_pieces], maxTime + 1)
times <- 1:T
deaths_length <- length(deaths)
deaths_starts <- seq(from = 1, to = deaths_length * 7, by = 7)
deaths_stops <- seq(from = 7, to = deaths_length * 7, by = 7)
twitter_symptons_length <- length(twitter_symptons)
stan_data <- list(
initial_time = initial_time,
n_beta_pieces = n_beta_pieces,
beta_left_t = beta_left_t,
beta_right_t = beta_right_t,
T = maxTime,
times = times,
n_disease_states = n_disease_states,
population = population,
deaths_length = length(deaths),
deaths_starts = deaths_starts,
deaths_stops = deaths_stops,
deaths = deaths,
real_data_length = length(beta_left_t) + length(beta_right_t) + 1,
real_data = c(beta_left_t, beta_right_t, population),
integer_data_length = 5,
integer_data = c(
maxTime,
length(beta_left_t),
length(beta_right_t),
length(beta_left_t),
n_disease_states
),
n_rho_twitter_symptons_pieces = n_rho_twitter_symptons_pieces,
rho_twitter_symptons_left_t = rho_twitter_symptons_left_t,
rho_twitter_symptons_right_t = rho_twitter_symptons_right_t,
twitter_symptons_start = twitter_symptons_start,
twitter_symptons = twitter_symptons,
twitter_symptons_length = twitter_symptons_length,
compute_likelihood = 1
)Running the Stan model, NOT WORKING getting NA values in beta_right_t.
library("cmdstanr")
#reminder to set working directory to a sibling to this file and the 'stan' dir below
model <- cmdstan_model("stan/deaths_twitter.stan")
#fit = model$optimize(data = stan_data) # doesn't work but params check out
set.seed(0)
initf <- function() {
n_beta_pieces <- stan_data$n_beta_pieces
initial_state_raw = c(runif(1, min = 0.99999, max = 1.0), runif(1, min = 0.0, max = 1.0))
beta_left = exp(runif(n_beta_pieces, min = -2, max = 0.5))
beta_right = exp(runif(n_beta_pieces, min = -2, max = 0.5))
dL = runif(1, min = 3.5, max = 4.0)
dI = runif(1, min = 2.2, max = 2.6)
dT = runif(1, min = 11.0, max = 13.0)
omega = 1/(1 + exp(-runif(1, min = -5, max = -3)))
list(
initial_state_raw = initial_state_raw,
beta_left = beta_left,
beta_right = beta_right,
dL = dL,
dI = dI,
dT = dT,
omega = omega
)
}
num_warmup <- 10
num_samples <- 10
num_chains <- 1
# model$sample(data=stan_data,
# seed=999,
# chains=num_chains,
# init = initf,
# iter_warmup=num_warmup,
# iter_sampling=num_samples,
# output_dir="output")