## ## Tutorial: Working with Nodal Attributes in Network Models ## Day 3 | Network Modeling for Epidemics ## library("EpiModel") # Clear environment rm(list = ls()) # Set seed set.seed(12345) # Network model ----------------------------------------------------------- # Initialize network num.g1 <- num.g2 <- 250 nw <- network_initialize(n = num.g1 + num.g2) nw # Set group attributes group <- rep(1:2, times = c(num.g1, num.g2)) nw <- set_vertex_attribute(nw, attrname = "group", value = group) nw get_vertex_attribute(nw, "group") # Degree distributions deg.dist.g1 <- c(0.40, 0.55, 0.04, 0.01) deg.dist.g2 <- c(0.54, 0.31, 0.10, 0.05) pois.dists <- c(dpois(0:2, lambda = 0.66), ppois(2, lambda = 0.66, lower.tail = FALSE)) # Plot degree par(mar = c(3,3,2,1), mgp = c(2,1,0), mfrow = c(1,1)) barplot(cbind(deg.dist.g1, deg.dist.g2, pois.dists), beside = TRUE, ylim = c(0, 0.6), col = 1:4) legend("topright", legend = paste0("deg", 0:3), pch = 15, col = 1:4, cex = 0.9, bg = "white") # Check degree check_degdist_bal(num.g1, num.g2, deg.dist.g1, deg.dist.g2) # Define the formation model formation <- ~edges + degree(0:1, by = "group") + nodematch("group") # Input the appropriate target statistics for each term target.stats <- c(165, 100, 137.5, 135, 77.5, 0) # Parameterize the dissolulation model coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 10) coef.diss # Fit the model est <- netest(nw, formation, target.stats, coef.diss) # Model diagnostics dx <- netdx(est, nsims = 10, nsteps = 500, ncores = 5, nwstats.formula = ~edges + degree(0:3, by = "group")) dx plot(dx) # Re-check degree check_degdist_bal(num.g1, num.g2, deg.dist.g1, deg.dist.g2) # Epidemic model ---------------------------------------------------------- # Parameterize model param <- param.net(inf.prob = 0.2, inf.prob.g2 = 0.2, rec.rate = 0.02, rec.rate.g2 = 0.02) # Initial conditions init <- init.net(i.num = 10, i.num.g2 = 10, r.num = 0, r.num.g2 = 0) # Control settings control <- control.net(type = "SIR", nsims = 5, nsteps = 500, ncores = 5) # Run the network model simulation with netsim sim <- netsim(est, param, init, control) sim # Plot the output # Plot 1 par(mfrow = c(1, 1)) plot(sim, main = "Disease State Sizes") par(mfrow = c(1,2)) plot(sim, y = c("i.num", "i.num.g2"), popfrac = TRUE, qnts = FALSE, ylim = c(0, 0.4), legend = TRUE) plot(sim, y = c("si.flow", "si.flow.g2"), qnts = FALSE, ylim = c(0, 2.5), legend = TRUE) # Plot 2 par(mfrow = c(1, 1), mar = c(0, 0, 0, 0)) plot(sim, type = "network", col.status = TRUE, at = 50, sims = "mean", shp.g2 = "square") # Add standardized incidence rates to sim sim <- mutate_epi(sim, ir.g1 = (si.flow / s.num) * 100 * 52, ir.g2 = (si.flow.g2 / s.num.g2) * 100 * 52) sim # Plot standardized incidence rates par(mar = c(3,3,2,1)) plot(sim, y = c("ir.g1", "ir.g2"), legend = TRUE, main = "Standardized Incidence Rates by Group") plot(sim, y = c("ir.g1", "ir.g2"), legend = TRUE, mean.smooth = FALSE, qnts = FALSE, mean.lwd = 1, main = "Standardized Incidence Rates by Group (Non-Smoothed)") # Convert model to a data frame for further analysis df <- as.data.frame(sim, out = "mean") head(df, 10) # Extract cumulative incidence sum(df$si.flow, na.rm = TRUE) sum(df$si.flow.g2, na.rm = TRUE) # Time step at max incidence and value which.max(df$ir.g1) df$ir.g1[which.max(df$ir.g1)] # Plot max inciidence plot(sim, y = "ir.g1", mean.smooth = FALSE, qnts = FALSE, mean.lwd = 1, main = "Time and Value of Peak Incidence in Females (Group 1)") points(which.max(df$ir.g1), df$ir.g1[which.max(df$ir.g1)], cex = 3) # Extract transmission matrix tm <- get_transmat(sim) head(tm, 15) mean(tm$infDur == 1) # Extract phylo trees tmPhylo <- as.phylo.transmat(tm) # Plot phylo trees par(mfrow = c(2, 2), mar = c(0, 0, 0, 0)) plot(tmPhylo[[1]], show.node.label = TRUE, root.edge = TRUE, cex = 0.5) plot(tmPhylo[[2]], show.node.label = TRUE, root.edge = TRUE, cex = 0.5) plot(tmPhylo[[3]], show.node.label = TRUE, root.edge = TRUE, cex = 0.5) plot(tmPhylo[[4]], show.node.label = TRUE, root.edge = TRUE, cex = 0.5) # Extract single tree names(tmPhylo) tmPhylo5 <- tmPhylo[[5]] # Plot single tree par(mfrow = c(2, 2), mar = c(0, 0, 0, 0)) plot(tmPhylo5, type = "phylogram", show.node.label = TRUE, root.edge = TRUE, cex = 0.5) plot(tmPhylo5, type = "cladogram", show.node.label = TRUE, root.edge = TRUE, cex = 0.5) plot(tmPhylo5, type = "fan", show.node.label = TRUE, root.edge = TRUE, cex = 0.5) plot(tmPhylo5, type = "unrooted", show.node.label = TRUE, root.edge = TRUE, cex = 0.5)