SEIRfansy
This R
package fits Extended Susceptible-Exposed-Infected-Recovery
(SEIR) Models for handling high false negative rate and symptom based
administration of diagnostic tests.
If the devtools package is not yet installed, install it first:
install.packages('devtools')
# install SEIRfansy from Github:
devtools::install_github('umich-biostatistics/SEIRfansy')
Once installed, load the package:
library(SEIRfansy)
#> Registered S3 methods overwritten by 'car':
#> method from
#> influence.merMod lme4
#> cooks.distance.influence.merMod lme4
#> dfbeta.influence.merMod lme4
#> dfbetas.influence.merMod lme4
For this example, we use the built-in package data set covid19
, which
contains dailies and totals of cases, recoveries, and deaths from the
COVID-19 outbreak in India from January 30 to September 21 of 2020.
You will need the dplyr
package for this example.
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
Training data set:
For training data, we use cases from April 1 to June 30
train = covid19[which(covid19$Date == "01 April "):which(covid19$Date == "30 June "),]
Testing data set:
For testing data, we use cases from July 1 to July 31
test = covid19[which(covid19$Date == "01 July "):which(covid19$Date == "31 July "),]
Data format for multinomial and Poisson distribution:
train_multinom =
train %>%
rename(Confirmed = Daily.Confirmed,
Recovered = Daily.Recovered,
Deceased = Daily.Deceased) %>%
dplyr::select(Confirmed, Recovered, Deceased)
test_multinom =
test %>%
rename(Confirmed = Daily.Confirmed,
Recovered = Daily.Recovered,
Deceased = Daily.Deceased) %>%
dplyr::select(Confirmed, Recovered, Deceased)
train_pois =
train %>%
rename(Confirmed = Daily.Confirmed) %>%
dplyr::select(Confirmed)
Initialize parameters:
N = 1341e6 # population size of India
data_initial = c(2059, 169, 58, 424, 9, 11)
pars_start = c(c(1,0.8,0.6,0.4,0.2), c(0.2,0.2,0.2,0.25,0.2))
phases = c(1,15,34,48,62)
If interest is in model estimation but not prediction, then use
SEIRfansy()
. Otherwise, use SEIRfansy.predict()
(see below).
?SEIRfansy
cov19est = SEIRfansy(data = train_multinom, init_pars = pars_start,
data_init = data_initial, niter = 1e3, BurnIn = 1e2,
model = "Multinomial", N = N, lambda = 1/(69.416 * 365),
mu = 1/(69.416 * 365), period_start = phases, opt_num = 1,
auto.initialize = TRUE, f = 0.15)
#> Finding MLE
#> 1 MLE run finished!
#>
#> MLE estimates :
#> beta = ( 0.18, 0.1, 0.25, 0.18, 0.18 )
#> r = ( 0.112, 0.531, 0.544, 0.65, 0.688 )
#>
#> MCMC:
#> Iter 100 A = 0 : 0.1772 0.1095 0.2517 0.1817 0.1749 0.1116 0.5091 0.5628
#> 0.6397 0.6707
#> Iter 200 A = 4e-04 : 0.1737 0.1133 0.2494 0.1829 0.1758 0.1117 0.4871 0.573
#> 0.6399 0.6686
#> Iter 300 A = 0.0026 : 0.1695 0.1169 0.2439 0.1847 0.177 0.1089 0.4676 0.5837
#> 0.6392 0.668
#> Iter 400 A = 9.7852 : 0.1669 0.1206 0.2406 0.189 0.1753 0.1101 0.4521 0.5945
#> 0.6368 0.6587
#> Iter 500 A = 0 : 0.1677 0.1241 0.232 0.1904 0.1777 0.1121 0.4379 0.597
#> 0.6393 0.652
#> Iter 600 A = 0 : 0.1688 0.1271 0.2261 0.19 0.1775 0.1147 0.4269 0.5992
#> 0.6358 0.6582
#> Iter 700 A = 0 : 0.1698 0.1299 0.2167 0.1931 0.1774 0.1163 0.4169 0.5985
#> 0.6415 0.6561
#> Iter 800 A = 96.4675 : 0.1682 0.134 0.2131 0.1941 0.1767 0.1155 0.4023
#> 0.6021 0.6432 0.6468
#> Iter 900 A = 0 : 0.1666 0.1409 0.2037 0.1955 0.178 0.1165 0.3833 0.5982
#> 0.6422 0.6415
#> Iter 1000 A = 0 : 0.1684 0.1432 0.1976 0.196 0.1792 0.1176 0.3724 0.5994
#> 0.6399 0.6398
#> Iter 1100 A = 4e-04 : 0.169 0.1467 0.1945 0.1969 0.1791 0.1154 0.3597 0.5875
#> 0.6358 0.6282
Inspect the results:
names(cov19est)
class(cov19est$mcmc_pars)
names(cov19est$plots)
Plot the results:
plot(cov19est, type = "trace")
plot(cov19est, type = "boxplot")
If interest is in model estimation and prediction, then use
SEIRfansy.predict()
, which first runs SEIRfansy()
internally, and
then predicts.
?SEIRfansy.predict
cov19pred = SEIRfansy.predict(data = train_multinom, init_pars = pars_start,
data_init = data_initial, T_predict = 60, niter = 1e3,
BurnIn = 1e2, data_test = test_multinom, model = "Multinomial",
N = N, lambda = 1/(69.416 * 365), mu = 1/(69.416 * 365),
period_start = phases, opt_num = 1,
auto.initialize = TRUE, f = 0.15)
#> Estimating ...
#>
#> Finding MLE
#> 1 MLE run finished!
#>
#> MLE estimates :
#> beta = ( 0.18, 0.1, 0.25, 0.18, 0.18 )
#> r = ( 0.112, 0.531, 0.544, 0.65, 0.688 )
#>
#> MCMC:
#> Iter 100 A = 0 : 0.176 0.1073 0.2544 0.1814 0.1747 0.11 0.5132 0.5471 0.6516
#> 0.6921
#> Iter 200 A = 1.683785e+12 : 0.1725 0.1105 0.2518 0.1851 0.1755 0.1104 0.4926
#> 0.5624 0.6498 0.6774
#> Iter 300 A = 0 : 0.1696 0.1153 0.2439 0.1879 0.1762 0.1119 0.4711 0.5755
#> 0.6533 0.6663
#> Iter 400 A = 0 : 0.1712 0.1188 0.2364 0.188 0.177 0.1119 0.4573 0.5834
#> 0.6465 0.6702
#> Iter 500 A = 0 : 0.1715 0.1215 0.2279 0.188 0.1772 0.111 0.4472 0.5995
#> 0.6534 0.6777
#> Iter 600 A = 498.6388 : 0.1699 0.1239 0.2246 0.1926 0.1756 0.1132 0.435
#> 0.6094 0.6533 0.6732
#> Iter 700 A = 0.0217 : 0.1689 0.1283 0.2174 0.1941 0.1781 0.1128 0.4204
#> 0.6056 0.655 0.6553
#> Iter 800 A = 407.4618 : 0.1703 0.1327 0.2114 0.1931 0.179 0.1147 0.4011
#> 0.6058 0.6508 0.6497
#> Iter 900 A = 1947.649 : 0.1703 0.1363 0.2038 0.195 0.1782 0.1166 0.3894
#> 0.6122 0.6487 0.6557
#> Iter 1000 A = 0 : 0.1696 0.1394 0.2005 0.1945 0.1784 0.1156 0.3812 0.61
#> 0.6448 0.6528
#> Iter 1100 A = 0 : 0.1737 0.1424 0.1936 0.1941 0.18 0.1164 0.3638 0.6017
#> 0.6499 0.6422
#>
#> Predicting ...
Inspect the results:
names(cov19pred)
class(cov19pred$prediction)
class(cov19pred$mcmc_pars)
names(cov19pred$plots)
Plot the results:
plot(cov19pred, type = "trace")
plot(cov19pred, type = "boxplot")
plot(cov19pred, type = "panel")
plot(cov19pred, type = "cases")
Ritwik Bhaduri, Ritoban Kundu, Soumik Purkayastha, Mike Kleinsasser, Lauren J Beesley, Bhramar Mukherjee. “EXTENDING THE SUSCEPTIBLE-EXPOSED-INFECTED-REMOVED(SEIR) MODEL TO HANDLE THE HIGH FALSE NEGATIVE RATE AND SYMPTOM-BASED ADMINISTRATION OF COVID-19 DIAGNOSTIC TESTS: SEIR-fansy.” medRxiv 2020.09.24.20200238; doi: https://doi.org/10.1101/2020.09.24.20200238
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.