load package
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.1, created on 2015-10-06
## 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'
Construct empty network
nw = network.initialize (n=150, directed =FALSE)
nw = set.vertex.attribute (nw, "race", rep(0:1, each = 75))
specify partnership
formation = ~edges + nodefactor("race") + nodematch("race") + concurrent
calculate target statistics
target.stats = c(30, 30, 23, 25)
construct dissolution model
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
model fit
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
construct network model
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 1.464
## Step length converged once. Increasing MCMC sample size.
## Iteration 2 of at most 20:
## The log-likelihood improved by 0.0729
## Step length converged twice. Stopping.
##
## This model was fit using MCMC. To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.
model diagnosis
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
## |*****|
##
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 35.183 0.173 9.256
## nodefactor.race.0 NA 32.020 NA 12.801
## nodefactor.race.1 30 38.347 0.278 10.220
## nodematch.race 23 26.847 0.167 7.777
## concurrent 25 28.330 0.133 8.492
##
## Dissolution Diagnostics
## -----------------------
## Target Sim Mean Pct Diff Sim SD
## Edge Duration 15.000 12.862 -0.143 12.413
## Pct Edges Diss 0.067 0.066 -0.011 0.043
Plot diagnosis#1
par(mar = c(3,3,1,1), mgp = c(2,1,0))
plot(dx)
Another plot diagnosis #2
par(mfrow = c(1, 2))
plot(dx, type = "duration")
plot(dx, type = "dissolution")