Interweaving vital and relational dynamics introduces multiple issues. Some of these have been known about for a while; others have been brought to light by the recent efforts to enable stergm to simulate off of the end of a networkDynamic object, given all of the latter's new functionality (e.g. net.obs.period). The presence of multiple similar and interrelated issues has at times made it hard to discuss and diagnose them; this is one attempt to tease them apart.
The size correction of Krivitsky et al. (2011) handles the issue of ergm parameters for changing population size; it does not deal with the additional issues discussed here
Steve's death correction handles (approximately, and well enough in practice so far) the issue of deaths causing relationships to break in ways not accounted for during estimation, which then causes durations to shorten and density to decline. This is regardless of the order in which events occur within a time step. We will call this the “Death Correction” for the rest of the document.
Because of the mechanics of both networkDynamic and stergm, Skye has suggested that the best way to have them interact (out of an imperfect set of options) is to always do relational dynamics first, then vital dynamics.
Unfortunately, it seems that #3 will require an additional correction beyond #1 and #2 above. This is because it introduces the additional issue that arrivals have to wait a time step to form ties. Neither of the corrections in #1 or #2 take care of this.
I have six examples below, designed to illustrate this. They include an ergm, a stergm without death correction, and a stergm with death correction. Each of these three is done with vital dynamics before relational dynamics, and then again vice versa. All models are Bernoulli.
A few things to keep in mind in them:
In order to isolate points #2 and #4 above, in all examples departures and arrivals are set equal to one another. This means that population size stays the same. To keep the code simple, I take account of this by simply selecting the nodes to die randomly, and removing all their ties. No need to actually remove them and then reintroduce an equal number of nodes, as this entails lots of bookkeeping but has the same effective outcome.
Given how the mechanics for simulating a stergm off of the end of a networkDynamic object have been in flux, and cause their own mental acrobatics, I have avoided this. Instead, for the stergm examples, I take a cross-sectional network object, simulate one time step from it, extract the network at the next time step as a new cross-sectional network object using network.collapse, rinse and repeat.
Instead of using the dissolution approximation (once called the Carnegie Correction), I use the exact value as presented in Carnegie et al (in press), which we know because this is a Bernouli model. This is seen below when I calculate stergm.coef.form. This means that we do not need to worry about that approximation being off and causing another layer of problems.
n <- 1000
mean.deg <- 2
duration <- 10
death.rate <- 0.05
time.steps <- 50
library(tergm)
mynet <- network.initialize(n, directed = F)
num.edges <- round(n * mean.deg/2)
mynet <- san(mynet ~ edges, target.stats = num.edges)
num.dying <- round(death.rate * n)
density <- num.edges/(n * (n - 1)/2)
ergm.coef <- ergm(mynet ~ edges)$coef
stergm.coef.diss <- log(duration - 1)
stergm.coef.form <- -log((1 + exp(stergm.coef.diss)) * (1 - density)/density +
1)
prob.diss.wo.death <- 1/duration
prob.neither.die <- (1 - death.rate)^2
prob.either.die <- 1 - prob.neither.die
logit <- function(x) log(x/(1 - x))
stergm.coef.diss.corrected <- logit(1 - (prob.diss.wo.death - prob.either.die)/prob.neither.die)
mynet <- san(mynet ~ edges, target.stats = num.edges)
meandeg.ergm.dt <- vector("numeric", time.steps)
meandeg.ergm.dt[1] <- summary(mynet ~ edges) * 2/n
for (i in 2:time.steps) {
dying <- sample(1:n, num.dying)
edges.to.rm <- unlist(lapply(dying, function(x) get.edgeIDs(mynet, x)))
mynet <- delete.edges(mynet, edges.to.rm)
mynet <- simulate(mynet ~ edges, coef = ergm.coef, control = control.simulate.formula(MCMC.burnin = 1e+05))
meandeg.ergm.dt[i] <- summary(mynet ~ edges) * 2/n
cat("Done ergm run ", i, "\n")
}
mynet <- san(mynet ~ edges, target.stats = num.edges)
meandeg.ergm.td <- vector("numeric", time.steps)
meandeg.ergm.td[1] <- summary(mynet ~ edges) * 2/n
for (i in 2:time.steps) {
mynet <- simulate(mynet ~ edges, coef = ergm.coef, control = control.simulate.formula(MCMC.burnin = 1e+05))
dying <- sample(1:n, num.dying)
edges.to.rm <- unlist(lapply(dying, function(x) get.edgeIDs(mynet, x)))
mynet <- delete.edges(mynet, edges.to.rm)
meandeg.ergm.td[i] <- summary(mynet ~ edges) * 2/n
cat("Done ergm run ", i, "\n")
}
mynet <- san(mynet ~ edges, target.stats = num.edges)
meandeg.stergm.dt.nc <- vector("numeric", time.steps)
meandeg.stergm.dt.nc[1] <- summary(mynet ~ edges) * 2/n
for (i in 2:time.steps) {
dying <- sample(1:n, num.dying)
edges.to.rm <- unlist(lapply(dying, function(x) get.edgeIDs(mynet, x)))
mynet <- delete.edges(mynet, edges.to.rm)
mynet <- simulate(mynet, formation = ~edges, dissolution = ~edges, coef.form = stergm.coef.form,
coef.diss = stergm.coef.diss, control = control.simulate.network(MCMC.burnin.min = 1e+05))
mynet <- network.collapse(mynet, at = 1)
meandeg.stergm.dt.nc[i] <- summary(mynet ~ edges) * 2/n
cat("Done ergm run ", i, "\n")
}
mynet <- san(mynet ~ edges, target.stats = num.edges)
meandeg.stergm.td.nc <- vector("numeric", time.steps)
meandeg.stergm.td.nc[1] <- summary(mynet ~ edges) * 2/n
for (i in 2:time.steps) {
mynet <- simulate(mynet, formation = ~edges, dissolution = ~edges, coef.form = stergm.coef.form,
coef.diss = stergm.coef.diss, control = control.simulate.network(MCMC.burnin.min = 1e+05))
dying <- sample(1:n, num.dying)
edges.to.rm <- unlist(lapply(dying, function(x) get.edgeIDs(mynet, x)))
mynet <- delete.edges(mynet, edges.to.rm)
mynet <- network.collapse(mynet, at = 1)
meandeg.stergm.td.nc[i] <- summary(mynet ~ edges) * 2/n
cat("Done ergm run ", i, "\n")
}
mynet <- san(mynet ~ edges, target.stats = num.edges)
meandeg.stergm.dt.c <- vector("numeric", time.steps)
meandeg.stergm.dt.c[1] <- summary(mynet ~ edges) * 2/n
for (i in 2:time.steps) {
dying <- sample(1:n, num.dying)
edges.to.rm <- unlist(lapply(dying, function(x) get.edgeIDs(mynet, x)))
mynet <- delete.edges(mynet, edges.to.rm)
mynet <- simulate(mynet, formation = ~edges, dissolution = ~edges, coef.form = stergm.coef.form,
coef.diss = stergm.coef.diss.corrected, control = control.simulate.network(MCMC.burnin.min = 1e+05))
mynet <- network.collapse(mynet, at = 1)
meandeg.stergm.dt.c[i] <- summary(mynet ~ edges) * 2/n
cat("Done ergm run ", i, "\n")
}
mynet <- san(mynet ~ edges, target.stats = num.edges)
meandeg.stergm.td.c <- vector("numeric", time.steps)
meandeg.stergm.td.c[1] <- summary(mynet ~ edges) * 2/n
for (i in 2:time.steps) {
mynet <- simulate(mynet, formation = ~edges, dissolution = ~edges, coef.form = stergm.coef.form,
coef.diss = stergm.coef.diss.corrected, control = control.simulate.network(MCMC.burnin.min = 1e+05))
dying <- sample(1:n, num.dying)
edges.to.rm <- unlist(lapply(dying, function(x) get.edgeIDs(mynet, x)))
mynet <- delete.edges(mynet, edges.to.rm)
mynet <- network.collapse(mynet, at = 1)
meandeg.stergm.td.c[i] <- summary(mynet ~ edges) * 2/n
cat("Done ergm run ", i, "\n")
}
plot(meandeg.ergm.dt, ylim = c(mean.deg * 0.25, mean.deg * 2), type = "l", lwd = 2,
xlab = "time step", ylab = "mean degree")
lines(c(1, time.steps), c(mean.deg, mean.deg), lwd = 2)
lines(meandeg.ergm.td, ylim = c(0, mean.deg * 2), lty = 2, lwd = 2)
lines(meandeg.stergm.dt.nc, ylim = c(0, mean.deg * 2), col = "blue", lwd = 2)
lines(meandeg.stergm.td.nc, ylim = c(0, mean.deg * 2), col = "blue", lty = 2,
lwd = 2)
lines(meandeg.stergm.dt.c, ylim = c(0, mean.deg * 2), col = "red", lwd = 2)
lines(meandeg.stergm.td.c, ylim = c(0, mean.deg * 2), col = "red", lty = 2,
lwd = 2)
legend(1, 4, c("ergm: Node then Edge Dynamics", "ergm: Edge then Node Dynamics",
"stergm w/o Death Correction: Node then Edge Dynamics", "stergm w/o Death Correction: Edge then Node Dynamics",
"stergm w/ Death Correction: Node then Edge Dynamics", "stergm w/ Death Correction: Edge then Node Dynamics"),
col = c("black", "black", "blue", "blue", "red", "red"), lty = c(1, 2, 1,
2, 1, 2), lwd = 3)