constraint not working in resimulation step

[ThomasKraft]

The MWE below demonstrates this issue. Given a hard constraint specified on ties within groups (vertex attribute), models behave as expected in producing only ties within groups. When a perturbation is added to the network attribute for groups in a duration = 1 model, my expectation is that during network resimulation ties would once again be restricted to the newly assigned groupings. Tracking this behavior in an example simulation however reveals that this does not work:

library(EpiModel)

num=300
toy <- network_initialize(n=num)
toy <- set_vertex_attribute(toy, "group", c(rep("1", num/3),rep("2", num/3),rep("3", num/3)))

ft <- netest(toy,
             formation = ~edges,
             constraints = ~blocks("group", levels2=~!diag(TRUE, length(unique(.nw%v%"group")))),
             target.stats = c(500),
             coef.diss = dissolution_coefs(dissolution = ~offset(edges), duration = 1), edapprox = T)

param <- param.net(inf.prob = 0.3, act.rate = 1)
init <- init.net(i.num = 10, r.num = 0)

move.net <- function(dat, at){
  net_summary <- summary(dat$nw[[1]]~edges + nodematch("group"), at = at)
  dat <- set_epi(dat, "proportion_before_move_nodematch", at, value=net_summary[2]/net_summary[1])

  current_groups <- get_attr(dat, "group")
  new_groups <- sample(current_groups)
  dat <- set_attr(dat, "group", new_groups)
  dat <- copy_datattr_to_nwattr(dat)

  net_summary <- summary(dat$nw[[1]]~edges + nodematch("group"), at = at)
  dat <- set_epi(dat, "proportion_after_move_nodematch", at, value=net_summary[2]/net_summary[1])
}

control <- control.net(type = NULL, nsteps = 50, nsims = 1, ncores = 1,
                       epi.by = "group",
                       move.FUN = move.net,
                       resim_nets.FUN = resim_nets,
                       infection.FUN = infection.net,
                       nwupdate.FUN = nwupdate.net, 
                       prevalence.FUN=prevalence.net,
                       tergmLite = F,
                       resimulate.network = T,
                       set.control.ergm=control.simulate.formula(MCMC.burnin=2e5, MCMC.interval=2e5),
                       module.order = c("move.FUN", "resim_nets.FUN", "infection.FUN", "nwupdate.FUN", "prevalence.FUN"))

sim <- netsim(ft, param, init, control)
sim$epi$proportion_before_move_nodematch  # would expect this to always be 1
sim$epi$proportion_after_move_nodematch   # would expect this to be <1

One possibility is that during the resimulation, this issue has to do with the use of the network at time at-1 as the time.start (but the basis seems like it should be fine). Is there an obvious way to achieve a result in which the constraint is upheld?

[chad-klumb]

The blocks constraint fixes the edge state of specific dyads: from ?InitErgmConstraint.blocks:

Any dyad whose toggle would produce a nonzero change statistic for a nodemix term with the same arguments will be fixed.

If you scramble the group attribute but leave the edges alone, you will introduce ties between groups. If the blocks constraint fixes the edge state of ties between groups, those ties will not dissolve in the next simulation step. If you want to scramble the group attribute and not have ties between groups, you should remove the between-group ties manually after scrambling the attribute.

[ThomasKraft]

Thanks @chad-klumb, this definitely solves it. Before closing this, can you just verify that the following is a reasonable approach for removing between group ties (section headed with "remove between-group ties")?

move.net <- function(dat, at){
  net_summary <- summary(dat$nw[[1]]~edges + nodematch("group"), at = at)
  dat <- set_epi(dat, "proportion_before_move_nodematch", at, value=net_summary[2]/net_summary[1])

  current_groups <- get_attr(dat, "group")
  new_groups <- current_groups
  new_groups[1:200] <- sample(current_groups[1:200])
  dat <- set_attr(dat, "group", new_groups)
  dat <- copy_datattr_to_nwattr(dat)

########### remove between-group ties
  moved <- which(new_groups != current_groups)
  ### deactivate those who moved
  dat$nw[[1]] <- deactivate.vertices(dat$nw[[1]], onset=at, terminus=Inf, deactivate.edges = TRUE, v = moved)
  ### reactivate those who moved
  dat$nw[[1]] <- activate.vertices(dat$nw[[1]], at=1, v = moved)
###########  

  net_summary <- summary(dat$nw[[1]]~edges + nodematch("group"), at = at)
  dat <- set_epi(dat, "proportion_after_move_nodematch", at, value=net_summary[2]/net_summary[1])
}

[chad-klumb]

I see the change to time.offset in 747a869d085f2ca0fbba7c2c166da9a01bb0ca1b makes this sort of thing awkward, which might be a good argument for using time.offset = 0 in netsim. I've put up the EpiModel@time_offset branch, which uses time.offset = 0 instead of time.offset = 1. With that branch, I would use the following code:

library(EpiModel)

num=300
toy <- network_initialize(n=num)
toy <- set_vertex_attribute(toy, "group", c(rep("1", num/3),rep("2", num/3),rep("3", num/3)))

ft <- netest(toy,
             formation = ~edges,
             constraints = ~blocks("group", levels2=~!diag(TRUE, length(unique(.nw%v%"group")))),
             target.stats = c(500),
             coef.diss = dissolution_coefs(dissolution = ~offset(edges), duration = 1), edapprox = T)

param <- param.net(inf.prob = 0.3, act.rate = 1)
init <- init.net(i.num = 10, r.num = 0)

move.net <- function(dat, at){
  net_summary <- summary(dat$nw[[1]]~edges + nodematch("group"), at = at)
  dat <- set_epi(dat, "proportion_before_move_nodematch", at, value=net_summary[2]/net_summary[1])

  current_groups <- get_attr(dat, "group")
  new_groups <- current_groups
  new_groups[1:200] <- sample(current_groups[1:200])
  dat <- set_attr(dat, "group", new_groups)
  dat <- copy_datattr_to_nwattr(dat)

########### remove between-group ties
  el <- as.edgelist(network.collapse(dat$nw[[1]], at = at))
  diffs <- which(new_groups[el[,1]] != new_groups[el[,2]])
  eld <- el[diffs,,drop=FALSE]
  eids <- unlist(get.dyads.eids(dat$nw[[1]], eld[,1], eld[,2]))
  if (length(eids) > 0) {
    dat$nw[[1]] <- deactivate.edges(dat$nw[[1]], onset = at, terminus = Inf, e = eids)
  }
###########  

  net_summary <- summary(dat$nw[[1]]~edges + nodematch("group"), at = at)
  dat <- set_epi(dat, "proportion_after_move_nodematch", at, value=net_summary[2]/net_summary[1])
}

control <- control.net(type = NULL, nsteps = 20, nsims = 1, ncores = 1,
                       epi.by = "group",
                       move.FUN = move.net,
                       resim_nets.FUN = resim_nets,
                       infection.FUN = infection.net,
                       nwupdate.FUN = nwupdate.net, 
                       prevalence.FUN=prevalence.net,
                       tergmLite = F,
                       resimulate.network = T,
                       set.control.ergm=control.simulate.formula(MCMC.burnin=2e5, MCMC.interval=2e5),
                       module.order = c("move.FUN", "resim_nets.FUN", "infection.FUN", "nwupdate.FUN", "prevalence.FUN"))

sim <- netsim(ft, param, init, control)
sim$epi$proportion_before_move_nodematch  # would expect this to always be 1
sim$epi$proportion_after_move_nodematch   # would expect this to always be 1

[ThomasKraft]

The EpiModel@time_offset branch is working great for this purpose. I will use that for now but would be curious how something like those changes could be integrated into a user option in the future (so that the functionality continues to be available with future EpiModel updates).

I'll also just note for posterity that in a model with demographic departures, adding retain.all.vertices=TRUE to the network.collapse() call above is needed.

[chad-klumb]

I think a single time scheme for networkDynamics would be preferable (changing it can affect the behavior of other modules). In light of this example I am inclining towards time.offset = 0, but the decision could bear some more thought.