| Title: | Predict Epidemic Curves with Generalized SEIR Modeling |
|---|---|
| Description: | Performs generalized Susceptible-Exposed-Infected-Recovered (SEIR) modeling to predict epidemic curves. The method is described in Peng et al. (2020) <doi:10.1101/2020.02.16.20023465>. |
| Authors: | Selcuk Korkmaz |
| Maintainer: | Selcuk Korkmaz <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.1.0 |
| Built: | 2026-10-05 09:53:11 UTC |
| Source: | https://github.com/selcukorkmaz/genSEIR |
This function compares the fitted and calcualted death and recovered ratios. The idea is to check whether the approximation of these ratios is appropriate.
checkRates(time, Q, R, D, kappaFun, lambdaFun, kappa, lambda, dt = 1/24)checkRates(time, Q, R, D, kappaFun, lambdaFun, kappa, lambda, dt = 1/24)
time |
time vector |
Q |
time histories of the quarantined/active cases |
R |
time histories of the recovered cases |
D |
time histories of the deceased cases |
kappaFun |
anonymous function approximating the death rate |
lambdaFun |
anonymous function approximating the recovery rate |
kappa |
mortality rate |
lambda |
cure rate |
dt |
a time step, default is 1/24. This oversample time to ensure that the algorithm converges. |
plots for death rate and recovery rate
Selcuk Korkmaz, [email protected]
Peng, L., Yang, W., Zhang, D., Zhuge, C., Hong, L. 2020. “Epidemic analysis of COVID-19 in China by dynamical modeling”, arXiv preprint arXiv:2002.06563.
Fit SEIQRDP function parameters used in the SEIQRDP function, used to model the time-evolution of an epidemic outbreak.
fit_SEIQRDP( Q, R, D, Npop, E0, I0, time, dt = 1/24, guess, ftol = 1e-06, ptol = 1e-06, gtol = 1e-06, epsfcn = 0.001, factor = 100, maxfev = 1000, maxiter = 100, nprint = 1, trace = TRUE, ... )fit_SEIQRDP( Q, R, D, Npop, E0, I0, time, dt = 1/24, guess, ftol = 1e-06, ptol = 1e-06, gtol = 1e-06, epsfcn = 0.001, factor = 100, maxfev = 1000, maxiter = 100, nprint = 1, trace = TRUE, ... )
Q |
time histories of the active cases |
R |
time histories of the recovered cases |
D |
time histories of the deceased cases |
Npop |
total population of the country |
E0 |
initial number of exposed cases |
I0 |
initial number of predicted infectious cases |
time |
a time vector |
dt |
the time step. This oversamples time to ensure that the algorithm converges |
guess |
initial guess parameters |
ftol |
nls.lm.control object. non-negative numeric. Default is |
ptol |
nls.lm.control object. non-negative numeric. Default is |
gtol |
nls.lm.control object. non-negative numeric. Default is |
epsfcn |
nls.lm.control object. Default is |
factor |
nls.lm.control object. Default is |
maxfev |
nls.lm.control object. Default is |
maxiter |
nls.lm.control object. Default is |
nprint |
nls.lm.control object. Default is |
trace |
set |
... |
further arguments |
a list of optimized parameters
Selcuk Korkmaz, [email protected]
Peng, L., Yang, W., Zhang, D., Zhuge, C., Hong, L. 2020. “Epidemic analysis of COVID-19 in China by dynamical modeling”, arXiv preprint arXiv:2002.06563.
start = "01/01/21" finish = "04/01/21" country = "Italy" dt = 1 f=30 covidData = getDataCOVID(start = start, finish = finish, country = country) Recovered = covidData$tableRecovered Deaths = covidData$tableDeaths Confirmed = covidData$tableConfirmed if(nrow(Recovered) == 1){ name = Recovered$CountryRegion }else{ name = paste0(Recovered$ProvinceState, " (",Recovered$CountryRegion,")") } recovered = Recovered[ ,5:ncol(covidData$tableRecovered)] deaths = Deaths[ ,5:ncol(covidData$tableDeaths)] confirmed = Confirmed[ ,5:ncol(covidData$tableConfirmed)] Npop = 60000000 alpha_guess = 0.05 beta_guess = 0.8 LT_guess = 7 Q_guess = 0.8 lambda_guess = c(0.01,0.001,10) kappa_guess = c(0.001,0.001,10) guess = c(alpha_guess, beta_guess, 1/LT_guess, Q_guess, lambda_guess, kappa_guess) Q0 = confirmed[1]-recovered[1]-deaths[1] I0 = 0.3*Q0 E0 = 0.3*Q0 R0 = recovered[1] D0 = deaths[1] Active = confirmed-recovered-deaths Active[Active<0] <- 0 Q=Active R=recovered D = deaths time = seq(as.Date(start, format = "%m/%d/%y"), as.Date(finish, format = "%m/%d/%y"), by = "1 day") params = fit_SEIQRDP(Q = Active, R = recovered, D = deaths, Npop = Npop, E0 = E0, I0 = I0, time = time, dt = dt, guess = guess, ftol = 1e-6, ptol = 1e-6, gtol = 1e-6, epsfcn = 0.001, factor = 100, maxfev = 1000,maxiter = 100, nprint = 1, trace = TRUE)start = "01/01/21" finish = "04/01/21" country = "Italy" dt = 1 f=30 covidData = getDataCOVID(start = start, finish = finish, country = country) Recovered = covidData$tableRecovered Deaths = covidData$tableDeaths Confirmed = covidData$tableConfirmed if(nrow(Recovered) == 1){ name = Recovered$CountryRegion }else{ name = paste0(Recovered$ProvinceState, " (",Recovered$CountryRegion,")") } recovered = Recovered[ ,5:ncol(covidData$tableRecovered)] deaths = Deaths[ ,5:ncol(covidData$tableDeaths)] confirmed = Confirmed[ ,5:ncol(covidData$tableConfirmed)] Npop = 60000000 alpha_guess = 0.05 beta_guess = 0.8 LT_guess = 7 Q_guess = 0.8 lambda_guess = c(0.01,0.001,10) kappa_guess = c(0.001,0.001,10) guess = c(alpha_guess, beta_guess, 1/LT_guess, Q_guess, lambda_guess, kappa_guess) Q0 = confirmed[1]-recovered[1]-deaths[1] I0 = 0.3*Q0 E0 = 0.3*Q0 R0 = recovered[1] D0 = deaths[1] Active = confirmed-recovered-deaths Active[Active<0] <- 0 Q=Active R=recovered D = deaths time = seq(as.Date(start, format = "%m/%d/%y"), as.Date(finish, format = "%m/%d/%y"), by = "1 day") params = fit_SEIQRDP(Q = Active, R = recovered, D = deaths, Npop = Npop, E0 = E0, I0 = I0, time = time, dt = dt, guess = guess, ftol = 1e-6, ptol = 1e-6, gtol = 1e-6, epsfcn = 0.001, factor = 100, maxfev = 1000,maxiter = 100, nprint = 1, trace = TRUE)
This function computes the matrix A that is found in: dY/dt = A*Y + F
getA(alpha, gamma, delta, lambda, kappa)getA(alpha, gamma, delta, lambda, kappa)
alpha |
protection rate |
gamma |
inverse of the average latent time |
delta |
rate of people entering in quarantine |
lambda |
cure rate |
kappa |
mortality rate |
The matrix A that is found in: dY/dt = A*Y + F
Selcuk Korkmaz, [email protected]
Peng, L., Yang, W., Zhang, D., Zhuge, C., Hong, L. 2020. “Epidemic analysis of COVID-19 in China by dynamical modeling”, arXiv preprint arXiv:2002.06563.
The function collects the updated COVID-19 data from the John Hopkins University.
getDataCOVID(country, start = NULL, finish = NULL)getDataCOVID(country, start = NULL, finish = NULL)
country |
name of the country. It should be a character string. |
start |
a start date in mm/dd/yy format. Start date can not be earlier than 01/22/20. Start date can not be later than finish date. If start date is |
finish |
a finish date in mm/dd/yy format. Finish date can not be earlier than start date. If finish date is |
a list of COVID-19 historical data including confirmed, death and recovered cases in desired time ranges.
Selcuk Korkmaz, [email protected]
Peng, L., Yang, W., Zhang, D., Zhuge, C., Hong, L. 2020. “Epidemic analysis of COVID-19 in China by dynamical modeling”, arXiv preprint arXiv:2002.06563.
covidData = getDataCOVID(country = "Italy", start = "05/01/20", finish = "12/31/20") recovered = covidData$tableRecovered deaths = covidData$tableDeaths confirmed = covidData$tableConfirmedcovidData = getDataCOVID(country = "Italy", start = "05/01/20", finish = "12/31/20") recovered = covidData$tableRecovered deaths = covidData$tableDeaths confirmed = covidData$tableConfirmed
This function provides a first estimate of the death rate, to faciliate convergence of the main algorithm.
getKappaFun( tTarget, Q, D, guess, ftol, ptol, gtol, epsfcn, factor, maxfev, maxiter, nprint, trace )getKappaFun( tTarget, Q, D, guess, ftol, ptol, gtol, epsfcn, factor, maxfev, maxiter, nprint, trace )
tTarget |
time vector |
Q |
target time-histories of the quarantined cases |
D |
target time-histories of the dead cases |
guess |
Initiail guess parameters for kappa |
ftol |
nls.lm.control object. non-negative numeric. Default is |
ptol |
nls.lm.control object. non-negative numeric. Default is |
gtol |
nls.lm.control object. non-negative numeric. Default is |
epsfcn |
nls.lm.control object. Default is |
factor |
nls.lm.control object. Default is |
maxfev |
nls.lm.control object. Default is |
maxiter |
nls.lm.control object. Default is |
nprint |
nls.lm.control object. Default is |
trace |
Set |
vector of estimation and optimization function for the death rate
Selcuk Korkmaz, [email protected]
Peng, L., Yang, W., Zhang, D., Zhuge, C., Hong, L. 2020. “Epidemic analysis of COVID-19 in China by dynamical modeling”, arXiv preprint arXiv:2002.06563.
This function provides a first estimate of the recovery rate, to faciliate convergence of the main algorithm.
getLambdaFun( tTarget, Q, R, guess, ftol, ptol, gtol, epsfcn, factor, maxfev, maxiter, nprint, trace )getLambdaFun( tTarget, Q, R, guess, ftol, ptol, gtol, epsfcn, factor, maxfev, maxiter, nprint, trace )
tTarget |
target time vector |
Q |
target time-histories of the quarantined cases |
R |
target time-histories of the recovered cases |
guess |
Initiail guess parameters for kappa |
ftol |
nls.lm.control object. non-negative numeric. Default is |
ptol |
nls.lm.control object. non-negative numeric. Default is |
gtol |
nls.lm.control object. non-negative numeric. Default is |
epsfcn |
nls.lm.control object. Default is |
factor |
nls.lm.control object. Default is |
maxfev |
nls.lm.control object. Default is |
maxiter |
nls.lm.control object. Default is |
nprint |
nls.lm.control object. Default is |
trace |
set |
vector of estimation and optimization function for the recovery rate
Selcuk Korkmaz, [email protected]
Peng, L., Yang, W., Zhang, D., Zhuge, C., Hong, L. 2020. “Epidemic analysis of COVID-19 in China by dynamical modeling”, arXiv preprint arXiv:2002.06563.
Anonymous function approximating the death rate
kappaFun(a, t)kappaFun(a, t)
a |
parameter vector |
t |
time vector |
No return value, called for side effects
Anonymous function approximating the recovery rate
lambdaFun(a, t)lambdaFun(a, t)
a |
parameter vector |
t |
time vector |
No return value, called for side effects
Model function
modelFun(Y, A, K)modelFun(Y, A, K)
Y |
time vector |
A |
the matrix A that is found in: dY/dt = A*Y + F |
K |
the zero matrix for the seven states |
No return value, called for side effects
This function creates plots for reported and predicted active, recovered and death cases.
plotSEIQRDP( object, active, recovered, deaths, title = NULL, params, checkRates = FALSE, ... )plotSEIQRDP( object, active, recovered, deaths, title = NULL, params, checkRates = FALSE, ... )
object |
a SEIQRDP result |
active |
time histories of the quarantined/active cases |
recovered |
time histories of the recovered cases |
deaths |
Time histories of the deceased cases |
title |
country name |
params |
fitted parameters by fit_SEIQRDP function |
checkRates |
if |
... |
other plot options |
plots for epidemic curves: active cases, recovered and deaths
Selcuk Korkmaz, [email protected]
Selcuk Korkmaz, [email protected]
start = "01/01/21" finish = "04/01/21" country = "Italy" dt = 1 f=30 covidData = getDataCOVID(start = start, finish = finish, country = country) Recovered = covidData$tableRecovered Deaths = covidData$tableDeaths Confirmed = covidData$tableConfirmed if(nrow(Recovered) == 1){ name = Recovered$CountryRegion }else{ name = paste0(Recovered$ProvinceState, " (",Recovered$CountryRegion,")") } recovered = Recovered[ ,5:ncol(covidData$tableRecovered)] deaths = Deaths[ ,5:ncol(covidData$tableDeaths)] confirmed = Confirmed[ ,5:ncol(covidData$tableConfirmed)] Npop = 60000000 alpha_guess = 0.05 beta_guess = 0.8 LT_guess = 7 Q_guess = 0.8 lambda_guess = c(0.01,0.001,10) kappa_guess = c(0.001,0.001,10) guess = c(alpha_guess, beta_guess, 1/LT_guess, Q_guess, lambda_guess, kappa_guess) Q0 = confirmed[1]-recovered[1]-deaths[1] I0 = 0.3*Q0 E0 = 0.3*Q0 R0 = recovered[1] D0 = deaths[1] Active = confirmed-recovered-deaths Active[Active<0] <- 0 Q=Active R=recovered D = deaths time = seq(as.Date(start, format = "%m/%d/%y"), as.Date(finish, format = "%m/%d/%y"), by = "1 day") params = fit_SEIQRDP(Q = Active, R = recovered, D = deaths, Npop = Npop, E0 = E0, I0 = I0, time = time, dt = dt, guess = guess, ftol = 1e-6, ptol = 1e-6, gtol = 1e-6, epsfcn = 0.001, factor = 100, maxfev = 1000,maxiter = 100, nprint = 1, trace = TRUE) res = SEIQRDP(alpha = params$alpha1, beta = params$beta1, gamma = params$gamma1, delta = params$delta1, lambda0 = params$Lambda1, kappa0 = params$Kappa1, Npop, E0, I0, Q0, R0, D0,lambdaFun = params$lambdaFun, kappaFun = params$kappaFun, tstart = start, tfinish = finish, dt = dt, f =f) p = plotSEIQRDP(res, Active, recovered, deaths, name, params, checkRates = TRUE) print(p)start = "01/01/21" finish = "04/01/21" country = "Italy" dt = 1 f=30 covidData = getDataCOVID(start = start, finish = finish, country = country) Recovered = covidData$tableRecovered Deaths = covidData$tableDeaths Confirmed = covidData$tableConfirmed if(nrow(Recovered) == 1){ name = Recovered$CountryRegion }else{ name = paste0(Recovered$ProvinceState, " (",Recovered$CountryRegion,")") } recovered = Recovered[ ,5:ncol(covidData$tableRecovered)] deaths = Deaths[ ,5:ncol(covidData$tableDeaths)] confirmed = Confirmed[ ,5:ncol(covidData$tableConfirmed)] Npop = 60000000 alpha_guess = 0.05 beta_guess = 0.8 LT_guess = 7 Q_guess = 0.8 lambda_guess = c(0.01,0.001,10) kappa_guess = c(0.001,0.001,10) guess = c(alpha_guess, beta_guess, 1/LT_guess, Q_guess, lambda_guess, kappa_guess) Q0 = confirmed[1]-recovered[1]-deaths[1] I0 = 0.3*Q0 E0 = 0.3*Q0 R0 = recovered[1] D0 = deaths[1] Active = confirmed-recovered-deaths Active[Active<0] <- 0 Q=Active R=recovered D = deaths time = seq(as.Date(start, format = "%m/%d/%y"), as.Date(finish, format = "%m/%d/%y"), by = "1 day") params = fit_SEIQRDP(Q = Active, R = recovered, D = deaths, Npop = Npop, E0 = E0, I0 = I0, time = time, dt = dt, guess = guess, ftol = 1e-6, ptol = 1e-6, gtol = 1e-6, epsfcn = 0.001, factor = 100, maxfev = 1000,maxiter = 100, nprint = 1, trace = TRUE) res = SEIQRDP(alpha = params$alpha1, beta = params$beta1, gamma = params$gamma1, delta = params$delta1, lambda0 = params$Lambda1, kappa0 = params$Kappa1, Npop, E0, I0, Q0, R0, D0,lambdaFun = params$lambdaFun, kappaFun = params$kappaFun, tstart = start, tfinish = finish, dt = dt, f =f) p = plotSEIQRDP(res, Active, recovered, deaths, name, params, checkRates = TRUE) print(p)
Runge-Kutta 4th Order Method to Solve Differential Equation
RK4(Y, A, K, dt)RK4(Y, A, K, dt)
Y |
initial values for seven states |
A |
the matrix A that is found in: dY/dt = A*Y + F |
K |
the zero matrix for the seven states |
dt |
the time step. This oversamples time to ensure that the algorithm converges |
ordinary differential equation result for the seven states
Selcuk Korkmaz, [email protected]
This function simulates the time-histories of an epidemic outbreak using a generalized SEIR model
SEIQRDP( alpha, beta, gamma, delta, lambda0, kappa0, Npop, E0, I0, Q0, R0, D0, lambdaFun, kappaFun, tstart, tfinish, dt = 1/24, f = 0 )SEIQRDP( alpha, beta, gamma, delta, lambda0, kappa0, Npop, E0, I0, Q0, R0, D0, lambdaFun, kappaFun, tstart, tfinish, dt = 1/24, f = 0 )
alpha |
fitted protection rate |
beta |
fitted infection rate |
gamma |
fitted Inverse of the average latent time |
delta |
fitted rate at which people enter in quarantine |
lambda0 |
fitted cure rate |
kappa0 |
fitted mortality rate |
Npop |
Total population of the sample |
E0 |
Initial number of exposed cases |
I0 |
Initial number of infectious cases |
Q0 |
Initial number of quarantined cases |
R0 |
Initial number of recovered cases |
D0 |
Initial number of dead cases |
lambdaFun |
anonymous function giving the time-dependant recovery rate |
kappaFun |
anonymous function giving the time-dependant death rate |
tstart |
start date |
tfinish |
finish date |
dt |
the time step. This oversamples time to ensure that the algorithm converges |
f |
future predictions |
a list of predicted cases including susceptible, exposed, infectious, quarantined, recovered, dead and insusceptible.
Selcuk Korkmaz, [email protected]
Selcuk Korkmaz, [email protected]
Peng, L., Yang, W., Zhang, D., Zhuge, C., Hong, L. 2020. “Epidemic analysis of COVID-19 in China by dynamical modeling”, arXiv preprint arXiv:2002.06563.
start = "01/01/21" finish = "04/01/21" country = "Italy" dt = 1 f=30 covidData = getDataCOVID(start = start, finish = finish, country = country) Recovered = covidData$tableRecovered Deaths = covidData$tableDeaths Confirmed = covidData$tableConfirmed if(nrow(Recovered) == 1){ name = Recovered$CountryRegion }else{ name = paste0(Recovered$ProvinceState, " (",Recovered$CountryRegion,")") } recovered = Recovered[ ,5:ncol(covidData$tableRecovered)] deaths = Deaths[ ,5:ncol(covidData$tableDeaths)] confirmed = Confirmed[ ,5:ncol(covidData$tableConfirmed)] Npop = 60000000 alpha_guess = 0.05 beta_guess = 0.8 LT_guess = 7 Q_guess = 0.8 lambda_guess = c(0.01,0.001,10) kappa_guess = c(0.001,0.001,10) guess = c(alpha_guess, beta_guess, 1/LT_guess, Q_guess, lambda_guess, kappa_guess) Q0 = confirmed[1]-recovered[1]-deaths[1] I0 = 0.3*Q0 E0 = 0.3*Q0 R0 = recovered[1] D0 = deaths[1] Active = confirmed-recovered-deaths Active[Active<0] <- 0 Q=Active R=recovered D = deaths time = seq(as.Date(start, format = "%m/%d/%y"), as.Date(finish, format = "%m/%d/%y"), by = "1 day") params = fit_SEIQRDP(Q = Active, R = recovered, D = deaths, Npop = Npop, E0 = E0, I0 = I0, time = time, dt = dt, guess = guess, ftol = 1e-6, ptol = 1e-6, gtol = 1e-6, epsfcn = 0.001, factor = 100, maxfev = 1000,maxiter = 100, nprint = 1, trace = TRUE) res = SEIQRDP(alpha = params$alpha1, beta = params$beta1, gamma = params$gamma1, delta = params$delta1, lambda0 = params$Lambda1, kappa0 = params$Kappa1, Npop, E0, I0, Q0, R0, D0,lambdaFun = params$lambdaFun, kappaFun = params$kappaFun, tstart = start, tfinish = finish, dt = dt, f =f)start = "01/01/21" finish = "04/01/21" country = "Italy" dt = 1 f=30 covidData = getDataCOVID(start = start, finish = finish, country = country) Recovered = covidData$tableRecovered Deaths = covidData$tableDeaths Confirmed = covidData$tableConfirmed if(nrow(Recovered) == 1){ name = Recovered$CountryRegion }else{ name = paste0(Recovered$ProvinceState, " (",Recovered$CountryRegion,")") } recovered = Recovered[ ,5:ncol(covidData$tableRecovered)] deaths = Deaths[ ,5:ncol(covidData$tableDeaths)] confirmed = Confirmed[ ,5:ncol(covidData$tableConfirmed)] Npop = 60000000 alpha_guess = 0.05 beta_guess = 0.8 LT_guess = 7 Q_guess = 0.8 lambda_guess = c(0.01,0.001,10) kappa_guess = c(0.001,0.001,10) guess = c(alpha_guess, beta_guess, 1/LT_guess, Q_guess, lambda_guess, kappa_guess) Q0 = confirmed[1]-recovered[1]-deaths[1] I0 = 0.3*Q0 E0 = 0.3*Q0 R0 = recovered[1] D0 = deaths[1] Active = confirmed-recovered-deaths Active[Active<0] <- 0 Q=Active R=recovered D = deaths time = seq(as.Date(start, format = "%m/%d/%y"), as.Date(finish, format = "%m/%d/%y"), by = "1 day") params = fit_SEIQRDP(Q = Active, R = recovered, D = deaths, Npop = Npop, E0 = E0, I0 = I0, time = time, dt = dt, guess = guess, ftol = 1e-6, ptol = 1e-6, gtol = 1e-6, epsfcn = 0.001, factor = 100, maxfev = 1000,maxiter = 100, nprint = 1, trace = TRUE) res = SEIQRDP(alpha = params$alpha1, beta = params$beta1, gamma = params$gamma1, delta = params$delta1, lambda0 = params$Lambda1, kappa0 = params$Kappa1, Npop, E0, I0, Q0, R0, D0,lambdaFun = params$lambdaFun, kappaFun = params$kappaFun, tstart = start, tfinish = finish, dt = dt, f =f)
Fitted Results for SEIQRDP
SEIQRDP_for_fitting(par, t, t0)SEIQRDP_for_fitting(par, t, t0)
par |
initial guess parameters |
t |
historical time vector |
t0 |
target time vector |
a data frame for fitted quarantined, recovered and deaths
Selcuk Korkmaz, [email protected]
Peng, L., Yang, W., Zhang, D., Zhuge, C., Hong, L. 2020. “Epidemic analysis of COVID-19 in China by dynamical modeling”, arXiv preprint arXiv:2002.06563.