Network Diffusion Modeling
WEN, Tzai-Hung (NTU Geography)
1. Basic ICMs with EpiModel: SI Model Stochasticity
stochastic individual contact models (ICMs) are designed to be agent-based microsimulation.
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/10mod## 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.flowsummary(mod, at = 125) # at time step 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. 258.3 114.450 0.516
## Infect. 242.7 114.450 0.484
## Total 501.0 0.000 1.000
## S -> I 5.0 2.667 NA
## ------------------------------head(as.data.frame(mod, out = "mean"))## time s.num i.num num si.flow
## 1 1 499.3 1.7 501 0.0
## 2 2 499.3 1.7 501 0.0
## 3 3 499.3 1.7 501 0.0
## 4 4 499.2 1.8 501 0.1
## 5 5 499.2 1.8 501 0.0
## 6 6 499.1 1.9 501 0.1tail(as.data.frame(mod, out = "vals", sim = 1))## time s.num i.num num si.flow
## 295 295 0 501 501 0
## 296 296 0 501 501 0
## 297 297 0 501 501 0
## 298 298 0 501 501 0
## 299 299 0 501 501 0
## 300 300 0 501 501 0Standard plotting output includes the means across simulations along with the inter-quartile range of values; both are smoothed values by default.
plot(mod)
plot(mod, y = "i.num", sim.lines = TRUE, mean.smooth = FALSE, qnts.smooth = FALSE)
2. Simulating Disease Epidemics over Dynamic Networks
Independent One-Mode SIS Model: a model of a Susceptible-Infected-Susceptible (SIS) epidemic in a closed population
2.1 Network Estimation and Diagnostics:
Specify a network structure, including features like size and compositional traits. Here, we construct an empty network of 1000 nodes with two races of equal node size.
nw <- network.initialize(n = 1000, directed = FALSE)
nw <- set.vertex.attribute(nw, "race", rep(0:1, each = 500))
- Model Parameters
formation <- ~edges + nodefactor("race") + nodematch("race") + concurrent
target.stats <- c(250, 375, 225, 100)
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 25)
coef.diss## Dissolution Coefficients
## =======================
## Dissolution Model: ~offset(edges)
## Edge Duration: 25
## Crude Coefficient: 3.178054
## Adjusted Coefficient: 3.178054
## Death rate: 0
- Model Fit
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 0.005906
## Step length converged once. Increasing MCMC sample size.
## Iteration 2 of at most 20:
## The log-likelihood improved by 0.008107
## Step length converged twice. Stopping.
##
## This model was fit using MCMC. To examine model diagnostics and check for degeneracy, use the mcmc.diagnostics() function.summary(est1)##
## ==========================
## Summary of model fit
## ==========================
##
## Formula: nw ~ edges + nodefactor("race") + nodematch("race") + concurrent
## <environment: 0x000000000d7d54b8>
##
## Iterations: 2 out of 20
##
## Monte Carlo MLE Results:
## Estimate Std. Error MCMC % p-value
## edges -9.8304138 0.2286235 0 <1e-04 ***
## nodefactor.race.1 0.6255021 0.0859362 0 <1e-04 ***
## nodematch.race 2.0113224 0.2196785 0 <1e-04 ***
## concurrent -0.0001165 0.1727395 0 0.999
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Log-likelihood was not estimated for this fit.
## To get deviances, AIC, and/or BIC from fit `object$fit` run
## > object$fit<-logLik(object$fit, add=TRUE)
## to add it to the object or rerun this function with eval.loglik=TRUE.
##
## Dissolution Coefficients
## =======================
## Dissolution Model: ~offset(edges)
## Edge Duration: 25
## Crude Coefficient: 3.178054
## Adjusted Coefficient: 3.178054
## Death rate: 0c.Model Diagnostics
dx <- netdx(est1, nsims = 5, nsteps = 500, 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: 500
##
## Formation Diagnostics
## -----------------------
## Target Sim Mean Pct Diff Sim SD
## edges 250 260.518 0.042 15.458
## nodefactor.race.0 NA 130.068 NA 15.011
## nodefactor.race.1 375 390.968 0.043 26.356
## nodematch.race 225 234.894 0.044 14.659
## concurrent 100 105.957 0.060 12.066
##
## Dissolution Diagnostics
## -----------------------
## Target Sim Mean Pct Diff Sim SD
## Edge Duration 25.00 23.646 -0.054 23.071
## Pct Edges Diss 0.04 0.040 0.003 0.012par(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")
2.2 Epidemic Simulation over Dynamic Networks
param <- param.net(inf.prob = 0.1, act.rate = 5, rec.rate = 0.02)
status.vector <- c(rbinom(500, 1, 0.1), rep(0, 500))
status.vector <- ifelse(status.vector == 1, "i", "s")
init <- init.net(status.vector = status.vector)
control <- control.net(type = "SIS", nsteps = 500, nsims = 10, epi.by = "race")
sim1 <- netsim(est1, param, init, control)Simulation output: Tables
sim1## EpiModel Simulation
## =======================
## Model class: netsim
##
## Simulation Summary
## -----------------------
## Model type: SIS
## No. simulations: 10
## No. time steps: 500
## No. NW modes: 1
##
## Model Parameters
## -----------------------
## inf.prob = 0.1
## act.rate = 5
## rec.rate = 0.02
##
## Model Output
## -----------------------
## Compartments: s.num s.num.race0 s.num.race1 i.num
## i.num.race0 i.num.race1 num num.race0 num.race1
## Flows: is.flow si.flow
## Networks: sim1 ... sim10
## Transmissions: sim1 ... sim10summary(sim1, at = 500)##
## EpiModel Summary
## =======================
## Model class: netsim
##
## Simulation Details
## -----------------------
## Model type: SIS
## No. simulations: 10
## No. time steps: 500
## No. NW modes: 1
##
## Model Statistics
## ------------------------------
## Time: 500
## ------------------------------
## mean sd pct
## Suscept. 498.4 27.633 0.498
## Infect. 501.6 27.633 0.502
## Total 1000.0 0.000 1.000
## S -> I 9.5 2.718 NA
## I -> S 9.2 2.098 NA
## ------------------------------head(as.data.frame(sim1), n=10)## time s.num s.num.race0 s.num.race1 i.num i.num.race0 i.num.race1 num
## 1 1 931.0 431.0 500.0 69.0 69.0 0.0 1000
## 2 2 925.7 427.1 498.6 74.3 72.9 1.4 1000
## 3 3 923.8 426.5 497.3 76.2 73.5 2.7 1000
## 4 4 921.8 425.3 496.5 78.2 74.7 3.5 1000
## 5 5 920.3 425.0 495.3 79.7 75.0 4.7 1000
## 6 6 919.5 425.1 494.4 80.5 74.9 5.6 1000
## 7 7 918.7 425.4 493.3 81.3 74.6 6.7 1000
## 8 8 918.9 426.0 492.9 81.1 74.0 7.1 1000
## 9 9 918.0 425.8 492.2 82.0 74.2 7.8 1000
## 10 10 917.6 426.8 490.8 82.4 73.2 9.2 1000
## num.race0 num.race1 is.flow si.flow
## 1 500 500 0.0 0.0
## 2 500 500 1.4 6.7
## 3 500 500 1.2 3.1
## 4 500 500 1.6 3.6
## 5 500 500 1.1 2.6
## 6 500 500 1.4 2.2
## 7 500 500 2.0 2.8
## 8 500 500 1.8 1.6
## 9 500 500 1.4 2.3
## 10 500 500 1.8 2.2head(as.data.frame(sim1, out = "vals", sim = 2)) # 2nd run## time s.num s.num.race0 s.num.race1 i.num i.num.race0 i.num.race1 num
## 1 1 931 431 500 69 69 0 1000
## 2 2 923 425 498 77 75 2 1000
## 3 3 921 424 497 79 76 3 1000
## 4 4 919 424 495 81 76 5 1000
## 5 5 919 425 494 81 75 6 1000
## 6 6 920 427 493 80 73 7 1000
## num.race0 num.race1 is.flow si.flow
## 1 500 500 0 0
## 2 500 500 1 9
## 3 500 500 0 2
## 4 500 500 0 2
## 5 500 500 1 1
## 6 500 500 4 3nw <- get_network(sim1, sim = 1)
nw## NetworkDynamic properties:
## distinct change times: 502
## maximal time range: -Inf until Inf
##
## Dynamic (TEA) attributes:
## Vertex TEAs: testatus.active
##
## Includes optional net.obs.period attribute:
## Network observation period info:
## Number of observation spells: 2
## Maximal time range observed: 1 until 501
## Temporal mode: discrete
## Time unit: step
## Suggested time increment: 1
##
## Network attributes:
## vertices = 1000
## directed = FALSE
## hyper = FALSE
## loops = FALSE
## multiple = FALSE
## bipartite = FALSE
## net.obs.period: (not shown)
## vertex.pid = vertex.pid
## edge.pid = edge.pid
## total edges= 5471
## missing edges= 0
## non-missing edges= 5471
##
## Vertex attribute names:
## active race testatus.active vertex.names vertex.pid
##
## Edge attribute names not shownSimulation output: Charts
plot(sim1)
plot(sim1, mean.line = FALSE, qnts = FALSE, sim.lines = TRUE)
plot(sim1, y = c("i.num.race0", "i.num.race1"), leg = TRUE)
Simulation output: Graphs
par(mfrow = c(1,2))
plot(sim1, type = "network", at = 1, col.status = TRUE, sims = 1, main = "Prevalence at t1")
plot(sim1, type = "network", at = 50, col.status = TRUE, sims = 1, main = "Prevalence at t50")
plot(sim1, type = "network", at = 100, col.status = TRUE, sims = 1, main = "Prevalence at t100")
plot(sim1, type = "network", at = 500, col.status = TRUE, sims = 1, main = "Prevalence at t500")
Export to images
setwd(“D:/R_Labs/PICs”)
m = 1
while (m <= 100){
filename = paste(“Epidemic_Day”, m, “.png”)
png(filename)
plot(sim1, type = “network”, at = m, col.status = TRUE, sims = 1)
title(paste(“Prevalence at Day”, m))
m = m + 1}