This short tutorial shows how to set birth or entry rates in a bipartite network model. Load EpiModel.

library(EpiModel)
if (packageVersion("EpiModel") != "1.2.5") {
  update.packages("EpiModel")
}

First, let’s set up and estimate a simple ERGM on a bipartite network, with equally sized modes (50 in mode 1 and 50 in mode 2). Using an edges term for the ERGM and a target statistic of 0.5 implies that there will be a mean degree of 1 for the network.

nw <- network.initialize(n = 100, bipartite = 50, directed = FALSE)
formation <- ~edges
target.stats <- 50
coef.diss <- dissolution_coefs(dissolution = ~offset(edges), duration = 20,
                               d.rate = 0.0021)
est <- netest(nw, formation, target.stats, coef.diss, verbose = FALSE)

In the first parameterization, we set the birth or entry rate of mode 1 with b.rate. The default setting for the birth or entry rate for mode 2 in bipartite models is NA. This means that the number of incoming nodes in mode 2 will not be a function of the size of mode 2, but of mode 1. We designed the default like this because we figured that people would use bipartite models to represent two-sex models, in which the first mode were the females and the second mode was the males. Similar to most demographic models for population growth, then the size of the incoming population of men is not based on the size of males but on the size of females due to the fertility process. This example shows the implications of this in an extreme form, in which we assume that the death or exit rate for mode 2 is twice that of females. Let’s run this and see what happens. Note you can lower the number of simulations for computational efficiency, but there is a lot of stochastic variation so we keep it at 16.

param <- param.net(inf.prob = 0.3, inf.prob.m2 = 0.15,
                   rec.rate = 0.02, rec.rate.m2 = 0.02,
                   b.rate = 0.01, b.rate.m2 = NA,
                   ds.rate = 0.01, ds.rate.m2 = 0.02,
                   di.rate = 0.01, di.rate.m2 = 0.02)
init <- init.net(i.num = 10, i.num.m2 = 10)
control <- control.net(type = "SIS", nsteps = 100, nsims = 16, ncores = 16)
mod1 <- netsim(est, param, init, control)
## 
## * Starting Network Simulation

Plotting the results, the left panel shows the size of the mode 1 and mode 2 populations over time in response to the nature birth and death processes. Unsurprisingly, the size of the mode 2 population declines over time as a reflection of the two-fold death rate. But interestingly, the birth rate for mode 2 equals that of mode 1 even when the mode 2 population size has considerably declined. Again, this is because the number of mode 2 births is a function only of the size of the first mode when the b.rate.m2 parameter is set to NA.

par(mfrow = c(1,2))
plot(mod1, y = c("num", "num.m2"), qnts = 0.5, 
     main = "Population Size", leg = TRUE)
plot(mod1, y = c("b.flow", "b.flow.m2"), qnts = 0.5, 
     ylim = c(0, 2), main = "Birth Rate", leg = TRUE)

Let’s compare that to an alternative model in which we explicitly set the birth rate of mode 2 to 0.01, the same as the birth rate of mode 1. What happens under the hood now? The birth rate of mode 2 will be based on the current size of mode 2, independent of mode 1. Now the “birth” process really reflects more of a general in-migration process where the migration rates of each mode are independent of one another. Under the same scenario where there is still a two-fold higher death or exit rate for mode 2, what will happen?

param <- param.net(inf.prob = 0.3, inf.prob.m2 = 0.15,
                   rec.rate = 0.02, rec.rate.m2 = 0.02,
                   b.rate = 0.01, b.rate.m2 = 0.01,
                   ds.rate = 0.01, ds.rate.m2 = 0.02,
                   di.rate = 0.01, di.rate.m2 = 0.02)
init <- init.net(i.num = 10, i.num.m2 = 10)
control <- control.net(type = "SIS", nsteps = 100, nsims = 16, ncores = 16)
mod2 <- netsim(est, param, init, control)
## 
## * Starting Network Simulation

This shows the expected phenomenon when the entry processes for the two modes are independent: as the population size of mode 2 declines, so does the “birth” rate.

par(mfrow = c(1,2))
plot(mod2, y = c("num", "num.m2"), qnts = 0.5, main = "Population Size", leg = TRUE)
plot(mod2, y = c("b.flow", "b.flow.m2"), qnts = 0.5, 
     ylim = c(0, 2), main = "Birth Rate", leg = TRUE)