data(concurrencyComparisonNets)
if (FALSE) {
# compare plots of each network at time 50
plot(network.extract(base,at=50),vertex.cex=0.5,edge.lwd=2)
plot(network.extract(monog,at=50),vertex.cex=0.5,edge.lwd=2)
plot(network.extract(middle,at=50),vertex.cex=0.5,edge.lwd=2)
# compare mean duration. These are at ~20, due to censoring
mean(as.data.frame(base)$duration)
mean(as.data.frame(middle)$duration)
mean(as.data.frame(monog)$duration)
# compare infection time series
plot(sapply(1:100,function(t){
sum(get.vertex.attribute.active(base,'status',at=t)==1)
}),col='black',xlab='time step', ylab='# infected'
)
points(sapply(1:100,function(t){
sum(get.vertex.attribute.active(monog,'status',at=t)==1)
}),col='blue')
points(sapply(1:100,function(t){
sum(get.vertex.attribute.active(middle,'status',at=t)==1)
}),col='red')
}
## The code below can be used generate the three example networks using EpiModel (v1.1.2)
## note that the networks have some attached simulation control variables deleted before
## being saved as the datasets. More recent versions of EpiModel use a different
## representation of the infection status variable.
if (FALSE) {
library(EpiModel)
# === example network parameters setup ===
params.base = list(
sim.length = 100,
size = 1000,
mean.deg = .75,
mean.rel.dur = 25,
net = network.initialize(1000, directed = F),
formation = ~edges,
dissolution = ~offset(edges)
)
params.middle = list(
sim.length = 100,
size = 1000,
mean.deg = .75,
mean.rel.dur = 25,
net = network.initialize(1000, directed = F),
formation = ~edges+concurrent,
dissolution = ~offset(edges),
targets = 90 # concurrent = 90
)
params.monog = list(
sim.length = 100,
size = 1000,
mean.deg = .75,
mean.rel.dur = 25,
net = network.initialize(1000, directed = F),
formation = ~edges+concurrent,
dissolution = ~offset(edges),
targets = 0 # concurrent = 0
)
# === function for estimating stergm, simulating network, and simulating epidemic ===
net.init <- function(params, nsims, adjust=F) {
for (name in names(params)) assign(name, params[[name]])
message('network init')
# generate initial network (defaults if not specified in params)
if (!exists('net', inherits=F)) {
net <- network.initialize(size,directed=F)
net
}
if (!exists('formation', inherits=F)) {
formation = ~edges
}
if (!exists('dissolution', inherits=F)) {
dissolution = ~offset(edges)
}
if (!is.null(mean.deg)) {
target.edges <- size/2 * mean.deg
density = target.edges / choose(size,2)
} else {
target.edges <- round(density*choose(size, 2))
}
print(target.edges)
# cludge to fix the monogamy bias in simulate
if (adjust) target.edges = target.edges*adjust
# target stats that does not include edges
if (exists('targets', inherits=F)) {
target.stats = c(target.edges, targets)
} else {
target.stats = target.edges
}
coef.diss <- dissolution_coefs(dissolution, mean.rel.dur)
# estimate the stergm
net.est = netest(net, formation, dissolution, target.stats, coef.diss
,set.control.ergm=control.ergm(MCMLE.maxit=200))
stats.form = update(formation, ~.+concurrent)
# simulate the dynamic network
#net.sim = netsim(net.est, nsteps = sim.length, nsims=nsims, stats.form=stats.form,
# control = control.simulate.network(MCMC.burnin.add=10))
# simulate the network dynamics and the epidemic
net.sim = netsim(net.est,
param.net(inf.prob=1),
init.net(i.num=1),
control.net(type='SI',
nsteps = sim.length,
nsims=nsims,
keep.network = TRUE)
)
#trans.sim = epiNet.simTrans(net.sim, "SI", vital=FALSE, i.num=1, trans.rate=1, tea=TRUE)
#print(summary(net.sim$stats[[1]]))
#plot(net.sim$stats[[1]][,'edges'], ylab='edges', xlab='time')
return(get_network(net.sim, sim = 1))
}
# === simulate example networks ===
base <- net.init(params.base, 1)
middle <- net.init(params.middle, 1)
monog <- net.init(params.monog, 1)
}
Run the code above in your browser using DataLab