Load add-on packages - deSolve - contains lsoda function - differential equation solver.
library (EpiModel)
## Loading required package: deSolve
##
## Attaching package: 'deSolve'
##
## The following object is masked from 'package:graphics':
##
## matplot
##
## Loading required package: networkDynamic
## Loading required package: network
## network: Classes for Relational Data
## Version 1.13.0 created on 2015-08-31.
## 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.8.0, created on 2015-09-14
## Copyright (c) 2015, 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: statnet.common
## Loading required package: ergm
##
## ergm: version 3.4.0, created on 2015-06-16
## Copyright (c) 2015, 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
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("ergm").
##
## NOTE: If you use custom ERGM terms based on 'ergm.userterms'
## version prior to 3.1, you will need to perform a one-time update
## of the package boilerplate files (the files that you did not write
## or modify) from 'ergm.userterms' 3.1 or later. See
## help('eut-upgrade') for instructions.
##
##
## Attaching package: 'ergm'
##
## The following objects are masked from 'package:network':
##
## as.edgelist, as.edgelist.matrix, as.edgelist.network
## Warning: replacing previous import by 'network::as.edgelist.network' when
## loading 'tergm'
## Warning: replacing previous import by 'network::as.edgelist.matrix' when
## loading 'tergm'
## Warning: replacing previous import by 'network::as.edgelist' when loading
## 'tergm'
##
## tergm: version 3.3.0, created on 2015-06-14
## Copyright (c) 2015, 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").
## Warning: replacing previous import by 'network::as.edgelist.network' when
## loading 'EpiModel'
## Warning: replacing previous import by 'network::as.edgelist.matrix' when
## loading 'EpiModel'
## Warning: replacing previous import by 'network::as.edgelist' when loading
## 'EpiModel'
param <- param.dcm(inf.prob = 0.2, act.rate = 0.25)
init <- init.dcm(s.num = 500, i.num = 1)
control <- control.dcm(type = "SI", nsteps = 500)
mod <- dcm(param, init, control)
mod
## EpiModel Simulation
## =======================
## Model class: dcm
##
## Simulation Summary
## -----------------------
## Model type: SI
## No. runs: 1
## No. time steps: 500
## No. groups: 1
##
## Model Parameters
## -----------------------
## inf.prob = 0.2
## act.rate = 0.25
##
## Model Output
## -----------------------
## Compartments: s.num i.num num
## Flows: si.flow
plot(mod)
summary(mod, at = 150)
## EpiModel Summary
## =======================
## Model class: dcm
##
## Simulation Summary
## -----------------------
## Model type: SI
## No. runs: 1
## No. time steps:
## No. groups: 1
##
## Model Statistics
## ------------------------------
## Time: 150 Run: 1
## ------------------------------
## n pct
## Suscept. 112.845 0.225
## Infect. 388.155 0.775
## Total 501.000 1.000
## S -> I 4.311 NA
## ------------------------------
param <- param.icm(inf.prob = 0.2, act.rate = 0.25)
init <- init.icm(s.num = 500, i.num = 1)
control <- control.icm(type = "SI", nsims = 10, nsteps = 300)
mod <- icm(param, init, control)
##
## * Starting ICM Simulation
## Sim = 1/10
## Sim = 2/10
## Sim = 3/10
## Sim = 4/10
## Sim = 5/10
## Sim = 6/10
## Sim = 7/10
## Sim = 8/10
## Sim = 9/10
## Sim = 10/10
mod
## EpiModel Simulation
## =======================
## Model class: icm
##
## Simulation Summary
## -----------------------
## Model type: SI
## No. simulations: 10
## No. time steps: 300
## No. groups: 1
##
## Model Parameters
## -----------------------
## inf.prob = 0.2
## act.rate = 0.25
##
## Model Output
## -----------------------
## Compartments: s.num i.num num
## Flows: si.flow
summary(mod, at = 125)
## EpiModel Summary
## =======================
## Model class: icm
##
## Simulation Details
## -----------------------
## Model type: SI
## No. simulations: 10
## No. time steps: 300
## No. groups: 1
##
## Model Statistics
## ------------------------------
## Time: 125
## ------------------------------
## mean sd pct
## Suscept. 289.3 63.540 0.577
## Infect. 211.7 63.540 0.423
## Total 501.0 0.000 1.000
## S -> I 5.6 2.875 NA
## ------------------------------
head(as.data.frame(mod, out = "mean"))
## time s.num i.num num si.flow
## 1 1 500.0 1.0 501 0.0
## 2 2 500.0 1.0 501 0.0
## 3 3 499.9 1.1 501 0.1
## 4 4 499.8 1.2 501 0.1
## 5 5 499.7 1.3 501 0.1
## 6 6 499.7 1.3 501 0.0
tail(as.data.frame(mod, out = "vals", sim = 1))
## time s.num i.num num si.flow
## 295 295 1 500 501 0
## 296 296 1 500 501 0
## 297 297 1 500 501 0
## 298 298 1 500 501 0
## 299 299 1 500 501 0
## 300 300 1 500 501 0
plot(mod)
plot(mod, y = "i.num", sim.lines = TRUE, mean.smooth = FALSE, qnts.smooth = FALSE)