diff --git a/2022-12-Multinets/README.md b/2022-12-Multinets/README.md new file mode 100644 index 0000000..da4407a --- /dev/null +++ b/2022-12-Multinets/README.md @@ -0,0 +1,37 @@ +# Epidemics with Multiple Networks + +## Description + +This example shows how to use core EpiModel functionality to include multiple, +interacting networks in a single epidemic model. The node set is required to be +the same for all networks. The edge formation and dissolution models may vary +from one network to the next, and the formation model for one network may depend +on the edge states in the other networks, represented via edge-dependent nodal +attributes. We give an example with two networks, including cross-network +dependency. The two network models are fit successively and then passed jointly +to `netsim`. + +### Modules + +We use built-in modules only, to illustrate their capabilities. + +### Parameters + +The cross-network dependency requires that we maintain the edge-dependent nodal +attributes during the `netsim` run via an appropriately specified `dat.updates` +argument to `control.net`. + +Additionally, certain arguments to `control.net` may vary from one network to +the next via the `multilayer` specification, and we illustrate this +functionality here by monitoring different statistics for the two networks. + +## Next Steps + +Multinetwork functionality is new to core EpiModel, and there are many possible +extensions. In the future it may be possible to have transmission probabilities +vary by edge type, for example. + +## Author + +Chad Klumb, University of Washington +(borrowing some code and comments from the SEIR example) diff --git a/2022-12-Multinets/model.R b/2022-12-Multinets/model.R new file mode 100644 index 0000000..e7ae97e --- /dev/null +++ b/2022-12-Multinets/model.R @@ -0,0 +1,107 @@ + +## +## Epidemics with Multiple Networks +## EpiModel Gallery (https://github.com/statnet/EpiModel-Gallery) +## +## Author: Chad Klumb (borrowing some code and comments from the SEIR example) +## Date: December 2022 +## + +## Load EpiModel +suppressMessages(library(EpiModel)) + +# Standard Gallery unit test lines +rm(list = ls()) +eval(parse(text = print(commandArgs(TRUE)[1]))) + +# Simulate data consistent with the target statistics we will use +nw <- network_initialize(500) +nw <- set_vertex_attribute(nw, "race", rep(c(0,1), length.out = 500)) +nw_san <- san(nw ~ edges + nodematch("race"), target.stats = c(90, 60)) +nw_san <- set_vertex_attribute(nw_san, "deg.net1", get_degree(nw_san)) +nw_san <- set_vertex_attribute(nw_san, "deg1+.net1", + pmin(get_degree(nw_san), 1)) +nw_san_2 <- san(nw_san ~ edges + degree(1) + nodefactor("deg1+.net1") + + Sum(matrix(seq_len(network.size(nw)), nrow = 1) ~ + degrange(1, by = "deg.net1", levels = -1),label="Sum"), + target.stats = c(75, 120, 10, 10)) +nw_san <- set_vertex_attribute(nw_san, "deg1+.net2", pmin(get_degree(nw_san_2), 1)) + +# confirm we have hit our intended targets for both edge types +summary(nw_san ~ edges + nodematch("race") + nodefactor("deg1+.net2")) +summary(nw_san_2 ~ edges + degree(1) + nodefactor("deg1+.net1")) + +# set edge-dependent nodal attributes on the network object for use in netest +nw <- set_vertex_attribute(nw, "deg1+.net1", pmin(get_degree(nw_san), 1)) +nw <- set_vertex_attribute(nw, "deg1+.net2", pmin(get_degree(nw_san_2), 1)) + +# Network model estimation ------------------------------------------------ + +# Define the formation models +formation.1 <- ~edges + nodematch("race") + nodefactor("deg1+.net2") +formation.2 <- ~edges + degree(1) + nodefactor("deg1+.net1") + +# Input the appropriate target statistics for each term +target.stats.1 <- c(90, 60, 10) +target.stats.2 <- c(75, 120, 10) + +# Parameterize the dissolution models +coef.diss.1 <- dissolution_coefs(dissolution = ~offset(edges), duration = 100) +coef.diss.2 <- dissolution_coefs(dissolution = ~offset(edges), duration = 75) + +# Fit the models +est.1 <- netest(nw, formation.1, target.stats.1, coef.diss.1) +est.2 <- netest(nw, formation.2, target.stats.2, coef.diss.2) + +# Model diagnostics +dx.1 <- netdx(est.1, nsims = 16, ncores = 4, nsteps = 500) +print(dx.1) +plot(dx.1) + +dx.2 <- netdx(est.2, nsims = 16, ncores = 4, nsteps = 500) +print(dx.2) +plot(dx.2) + +# Epidemic model simulation ----------------------------------------------- + +# Model parameters +param <- param.net(inf.prob = 0.5, act.rate = 2, + ei.rate = 0.01, ir.rate = 0.01) + +# Initial conditions +init <- init.net(i.num = 10) + +# Control settings +control <- control.net(type = "SI", + nsteps = 500, + nsims = 16, + ncores = 4, + tergmLite = TRUE, + resimulate.network = TRUE, + dat.updates = function(dat, at, network) { + if (network == 0) { + # update deg1+.net2 prior to network 1 simulation, + # as network 1's formation model depends on this + # nodal attribute + dat <- set_attr(dat, "deg1+.net2", pmin(get_degree(dat$el[[2]]), 1)) + } else if (network == 1) { + # update deg1+.net1 prior to network 2 simulation, + # as network 2's formation model depends on this + # nodal attribute + dat <- set_attr(dat, "deg1+.net1", pmin(get_degree(dat$el[[1]]), 1)) + } else if (network == 2) { + # update deg1+.net2 after network 2 simulation, + # for summary statistics + dat <- set_attr(dat, "deg1+.net2", pmin(get_degree(dat$el[[2]]), 1)) + } + return(dat) + }) + +# Run the network model simulation with netsim +sim <- netsim(list(est.1, est.2), param, init, control) + +# Examine results +print(sim, network = 1) +print(sim, network = 2) +plot(sim, network = 1, type = "formation") +plot(sim, network = 2, type = "formation")