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'

Set Network Structure

nw <- network.initialize(n = 150, directed = FALSE)
nw <- set.vertex.attribute(nw, "race", rep(0:1, each = 75))

Formulate Partnership Formation

formation <- ~edges + nodefactor("race") + nodematch("race") + concurrent

Calculate the target statistics

target.stats <- c(30, 30, 23, 25)
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 15)
coef.diss
## Dissolution Coefficients
## =======================
## Dissolution Model: ~offset(edges)
## Edge Duration: 15
## Crude Coefficient: 2.639057
## Adjusted Coefficient: 2.639057
## Death rate: 0

Fit Model

function (nw, formation, target.stats, coef.diss, constraints, 
    coef.form = NULL, edapprox = TRUE, output = "fit", set.control.ergm, 
    set.control.stergm, nonconv.error = FALSE, verbose = FALSE) 
NULL 
## function (nw, formation, target.stats, coef.diss, constraints, 
##     coef.form = NULL, edapprox = TRUE, output = "fit", set.control.ergm, 
##     set.control.stergm, nonconv.error = FALSE, verbose = FALSE) 
## NULL
est1 <- netest(nw, formation, target.stats, coef.diss, edapprox = TRUE)
## Starting maximum likelihood estimation via MCMLE:
## Iteration 1 of at most 20: 
## The log-likelihood improved by 2.689 
## Step length converged once. Increasing MCMC sample size.
## Iteration 2 of at most 20: 
## The log-likelihood improved by 0.6713 
## Step length converged twice. Stopping.
## 
## This model was fit using MCMC.  To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.

Simulate Model Fit Dynamics

dx <- netdx(est1, nsims = 5, nsteps = 120,
            nwstats.formula = ~edges + nodefactor("race", base = 0) +
                               nodematch("race") + concurrent)
## 
## Network Diagnostics
## -----------------------
## - Simulating 5 networks
##   |*****|
## - Calculating formation statistics
## - Calculating duration statistics
##   |*****|
## - Calculating dissolution statistics
##   |*****|
##

Print output

dx
## EpiModel Network Diagnostics
## =======================
## Diagnostic Method: Dynamic
## Simulations: 5
## Time Steps per Sim: 120
## 
## Formation Diagnostics
## ----------------------- 
##                   Target Sim Mean Pct Diff Sim SD
## edges                 30   36.127    0.204  6.682
## nodefactor.race.0     NA   42.850       NA  7.766
## nodefactor.race.1     30   29.403   -0.020 10.294
## nodematch.race        23   27.643    0.202  5.128
## concurrent            25   28.990    0.160  6.388
## 
## Dissolution Diagnostics
## ----------------------- 
##                Target Sim Mean Pct Diff Sim SD
## Edge Duration  15.000   12.881   -0.141 11.962
## Pct Edges Diss  0.067    0.067    0.007  0.044

Plot Diagnostics

par(mar = c(3,3,1,1), mgp = c(2,1,0))
plot(dx)

par(mfrow = c(1, 2))
plot(dx, type = "duration")
plot(dx, type = "dissolution")

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

Print

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.  250.8  126.453  0.501
## Infect.   250.2  126.453  0.499
## Total     501.0    0.000  1.000
## S -> I      4.3    2.263     NA
## ------------------------------
head(as.data.frame(mod, out = "mean"))
##   time s.num i.num num si.flow
## 1    1 499.5   1.5 501     0.0
## 2    2 499.5   1.5 501     0.0
## 3    3 499.5   1.5 501     0.0
## 4    4 499.5   1.5 501     0.0
## 5    5 499.5   1.5 501     0.0
## 6    6 499.3   1.7 501     0.2
head(as.data.frame(mod, out = "vals", sim = 1))
##   time s.num i.num num si.flow
## 1    1   499     2 501       0
## 2    2   499     2 501       0
## 3    3   499     2 501       0
## 4    4   499     2 501       0
## 5    5   499     2 501       0
## 6    6   499     2 501       0
plot(mod)

par(mfrow = c(1,1), mar = c(3,3,1,1), mgp = c(2,1,0))
plot(mod)

plot(mod, y = "i.num", sim.lines = TRUE, mean.smooth = FALSE, qnts.smooth = FALSE)

Plot Static Network of Different Time Steps of Different Simualtions

par(mfrow = c(1,2), mar = c(0,0,1,0))
plot(mod, type = "network", at = 1, col.status = TRUE,
     main = "Prevalence at t1")
plot(mod, type = "network", at = 120, col.status = TRUE,
     main = "Prevalence at t120")