SEIQRDP: Simulate generalized SEIR model

Description

This function simulates the time-histories of an epidemic outbreak using a generalized SEIR model

Usage

SEIQRDP(
  alpha,
  beta,
  gamma,
  delta,
  lambda0,
  kappa0,
  Npop,
  E0,
  I0,
  Q0,
  R0,
  D0,
  lambdaFun,
  kappaFun,
  tstart,
  tfinish,
  dt = 1/24,
  f = 0
)

Arguments

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

Initial number of exposed cases

Initial number of infectious cases

Initial number of quarantined cases

Initial number of recovered cases

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

the time step. This oversamples time to ensure that the algorithm converges

number of days for future predictions

Value

a list of predicted cases including susceptible, exposed, infectious, quarantined, recovered, dead and insusceptible.

References

Peng, L., Yang, W., Zhang, D., Zhuge, C., Hong, L. 2020. <U+201C>Epidemic analysis of COVID-19 in China by dynamical modeling<U+201D>, arXiv preprint arXiv:2002.06563.

https://www.mathworks.com/matlabcentral/fileexchange/74545-generalized-seir-epidemic-model-fitting-and-computation

Examples

Run this code

# NOT RUN {
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 = list(alpha_guess,
            beta_guess,
            1/LT_guess,
            Q_guess,
            lambda_guess[1],
            lambda_guess[2],
            lambda_guess[3],
            kappa_guess[1],
            kappa_guess[2],
            kappa_guess[3])

 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 = c(params$lambda01, params$lambda02, params$lambda03),
               kappa0 = c(params$kappa01, params$kappa02, params$kappa03),
               Npop, E0, I0, Q0, R0, D0,lambdaFun = params$lambdaFun,
               kappaFun = params$kappaFun, tstart = start, tfinish = finish,
               dt = dt, f =f)
# }