The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

R package SEIRfansy

Extended Susceptible-Exposed-Infected-Recovery Model

Overview

This R package fits Extended Susceptible-Exposed-Infected-Recovery (SEIR) Models for handling high false negative rate and symptom based administration of diagnostic tests.

Installation

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

Example Usage

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.

Setup

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)

SEIRfansy()

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")

SEIRfansy.predict()

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")

Current Suggested Citation

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

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.