EpiModel

EpiModel is an R package that provides tools for simulating and analyzing mathematical models of infectious disease dynamics. Supported epidemic model classes include

Disease types include SI, SIR, and SIS epidemics with and without demography, with utilities available for expansion to construct and simulate epidemic models of arbitrary complexity. The network model class is based on the statistical framework of temporal exponential random graph models (ERGMs) implementated in the Statnet suite of software for R.

Similar package to EpiModel is EpiDynamics.

Why EpiModel?

We can solve a simple SIR model with deSolve package in R. To do that we need

require(deSolve)
## Loading required package: deSolve
#time specification and state variables
times = seq(0, 26, by = 1/10)
parameters = c(mu = 0, N = 1, beta = 2, gamma = 1/2)
state = c(S = 0.999, I = 0.001, R = 0)

#model equations
sirmod = function (t, state, parameters){
  with(as.list(c(state, parameters)),{
    # model equations
    dS = mu*(N-S) - beta *S * I/N
    dI = beta *S * I/N - (mu+gamma)*I
    dR = gamma*I - mu*R
    result = c(dS, dI, dR)
    #return gradient
    list(result)
  })# end of with (as.list...
}

#Model integration
out=ode(y=state, times=times, func=sirmod, parms=parameters)
out=as.data.frame(out)

#Graphs
plot(x=out$time, y=out$S, ylab="Fraction", xlab=
"Time", type="l", col ="blue", lwd = 3)
lines(x=out$time, y=out$I, col="red",lwd = 3)
lines(x=out$time, y=out$R, col="green",lwd = 3)

#install.packages("EpiModel", dependencies = TRUE)
require(EpiModel)
## Loading required package: EpiModel
## Loading required package: networkDynamic
## Loading required package: network
## network: Classes for Relational Data
## Version 1.15 created on 2019-04-01.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
##                     Mark S. Handcock, University of California -- Los Angeles
##                     David R. Hunter, Penn State University
##                     Martina Morris, University of Washington
##                     Skye Bender-deMoll, University of Washington
##  For citation information, type citation("network").
##  Type help("network-package") to get started.
## 
## networkDynamic: version 0.10.0, created on 2019-04-04
## Copyright (c) 2019, Carter T. Butts, University of California -- Irvine
##                     Ayn Leslie-Cook, University of Washington
##                     Pavel N. Krivitsky, University of Wollongong
##                     Skye Bender-deMoll, University of Washington
##                     with contributions from
##                     Zack Almquist, University of California -- Irvine
##                     David R. Hunter, Penn State University
##                     Li Wang
##                     Kirk Li, University of Washington
##                     Steven M. Goodreau, University of Washington
##                     Jeffrey Horner
##                     Martina Morris, University of Washington
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("networkDynamic").
## Loading required package: tergm
## Loading required package: ergm
## 
## ergm: version 3.10.4, created on 2019-06-10
## Copyright (c) 2019, Mark S. Handcock, University of California -- Los Angeles
##                     David R. Hunter, Penn State University
##                     Carter T. Butts, University of California -- Irvine
##                     Steven M. Goodreau, University of Washington
##                     Pavel N. Krivitsky, University of Wollongong
##                     Martina Morris, University of Washington
##                     with contributions from
##                     Li Wang
##                     Kirk Li, University of Washington
##                     Skye Bender-deMoll, University of Washington
##                     Chad Klumb
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("ergm").
## NOTE: Versions before 3.6.1 had a bug in the implementation of the
## bd() constriant which distorted the sampled distribution somewhat.
## In addition, Sampson's Monks datasets had mislabeled vertices. See
## the NEWS and the documentation for more details.
## NOTE: Some common term arguments pertaining to vertex attribute
## and level selection have changed in 3.10.0. See terms help for
## more details. Use 'options(ergm.term=list(version="3.9.4"))' to
## use old behavior.
## 
## tergm: version 3.6.1, created on 2019-06-12
## Copyright (c) 2019, Pavel N. Krivitsky, University of Wollongong
##                     Mark S. Handcock, University of California -- Los Angeles
##                     with contributions from
##                     David R. Hunter, Penn State University
##                     Steven M. Goodreau, University of Washington
##                     Martina Morris, University of Washington
##                     Nicole Bohme Carnegie, New York University
##                     Carter T. Butts, University of California -- Irvine
##                     Ayn Leslie-Cook, University of Washington
##                     Skye Bender-deMoll
##                     Li Wang
##                     Kirk Li, University of Washington
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("tergm").
#shiny packages
#epiweb("dcm")

#a bit digression
  #epiweb("icm")
  
#epiweb("net")
#model parameters
param <- param.dcm(inf.prob = 0.2, act.rate = 1, rec.rate = 1/20,
                   a.rate = 1/95, ds.rate = 1/100, di.rate = 1/80,
                   dr.rate = 1/100)
#state variables - initial state
init <- init.dcm(s.num = 1000, i.num = 1, r.num = 0)
#model equations and time specification
control <- control.dcm(type = "SIR", nsteps = 500, dt = 0.5)
#model integration
mod <- dcm(param, init, control)

#Plotting
par(mar = c(3.2, 3, 2, 1), mgp = c(2, 1, 0), mfrow = c(1, 2))
plot(mod, popfrac = FALSE, alpha = 0.5,
     lwd = 4, main = "Compartment Sizes")
plot(mod, y = "si.flow", lwd = 4, col = "firebrick",
     main = "Disease Incidence", legend = "n")

#schematics
par(mfrow = c(1, 1))
comp_plot(mod, at = 50, digits = 1)

Extending to user’s dcm (SEIR)

#model equation
SEIRme <- function (t,t0, parms){
  with(as.list(c(t0,parms)),{
    #population size
    
    num = s.num +e.num + i.num + r.num
    
    #force of infection
    lambda <- inf.prob*act.rate* i.num/num 
    
    #SEIR model
    dS <- d.rate*(num*(1-p.vac)-s.num) - lambda*s.num
    dE <- lambda*s.num - (l.rate + d.rate)*e.num
    dI <- l.rate*e.num - (r.rate + di.rate + d.rate)*i.num
    dR <- r.rate*i.num - d.rate*r.num + d.rate*num*p.vac
  
  
  #compartments and flows are part of the derivatives vector
  list(c(dS,dE,dI,dR,
         se.flow = lambda*s.num,
         ei.flow = l.rate*e.num,
         ir.flow = r.rate*r.num,
         d.flow = di.rate*i.num),
         num = num,
         i.prev = i.num/num,
         ei.prev = (e.num+i.num)/num)
  })
}
#model parameters
param <- param.dcm(inf.prob = 0.2, act.rate = 1, r.rate = 1/20,
                   l.rate = 1/20,p.vac = 0.1, d.rate = 1/100,
                   di.rate = c(0.005,0.0125,0.01,0.05))
#state variables - initial state
init <- init.dcm(s.num = 1e6, e.num = 10, i.num = 0, r.num = 0,
                  se.flow = 0, ei.flow = 0, ir.flow = 0, d.flow = 0)
#model equations and time specification
control <- control.dcm(nsteps = 500, dt = 1, new.mod = SEIRme)
#model integration
mod <- dcm(param, init, control)
# 
# #Plotting
par(mar = c(3.2, 3, 2, 1), mgp = c(2, 1, 0), mfrow = c(1, 2))
plot(mod,y="i.num", popfrac = FALSE, alpha = 0.5,
     lwd = 4, main = "Prevalence")
plot(mod, y = "se.flow", lwd = 4,
     main = "Disease Incidence", legend = "n")

#

Interesting?

This is the epiModel official website:

http://www.epimodel.org/tut.html