# Module 10 Lab: Network Model Parameterization (Worked Solution) # # Accompanies the SISMID course "Network Modeling for Epidemics" (NME). # Source page: mod10-Lab.qmd # # Practice with the two ideas from the Module 10 slides: BALANCE (the number # of ties group A has with group B must equal the number B has with A) and # DEGREES OF FREEDOM (you can only specify as many target statistics as the # terms already in the model leave room for). We start from an egocentric # sample of 20 heterosexual respondents in two communities, reason our way # term by term to a formation model and its target statistics, fit and # diagnose the TERGM, then run an SIS epidemic on the simulated network. # # Two standing caveats. First, no one would ever fit a model from 20 # respondents; the sample is this small only so the arithmetic stays visible # by hand. Second, the ergm.ego package automates most of what we do below. # We do it manually here so the logic is explicit: so you can explain it to a # reviewer, and so you can calculate the statistics yourself when you need to. # # This is the worked solution. Try the specification yourself before reading # it; several of the decisions below have no single correct answer, and the # value of the exercise is in confronting them. # # Run top to bottom. set.seed(1) matches the seed on the course page, but # netdx() and netsim() below simulate in parallel, so your numbers will differ # slightly from the rendered page regardless of the seed. Set ncores to the # number of cores you actually have. # ---- Setup ---- # (We leave a workspace-clear commented out so that sourcing this script will # not erase your other objects. Uncomment it for a clean slate.) # rm(list = ls()) library(EpiModel) set.seed(1) # ===================================================================== # THE DATA # ===================================================================== # A sample of 20 heterosexuals living in two communities. Each was asked for # their sex and community, their number of current partners (all opposite # sex), and the community of each of those partners. NA in p1.community or # p2.community means the respondent did not report a partner in that slot. mydf <- data.frame( sex = c('F', 'M', 'F', 'M', 'F', 'M', 'M', 'F', 'F', 'F', 'F', 'M', 'F', 'F', 'M', 'F', 'M', 'M', 'M', 'M'), community = c(2,1,2,2,1,2,2,2,2,2, 1,1,1,1,2,2,1,1,2,2), num.partners = c(0,0,0,0,2,0,1,1,1,1, 0,2,1,1,2,1,1,2,0,1), p1.community = c(NA, NA, NA, NA, 1, NA, 2, 1, 2, 2, NA, 2, 1, 1, 2, 1, 2, 1, NA, 2), p2.community = c(NA, NA, NA, NA, 1, NA, NA, NA, NA, NA, NA, 2, NA, NA, 1, NA, NA, 1, NA, NA) ) mydf # The sample was constructed to be convenient: equal numbers of females and # males, with the same community split within each sex. That is, 20% community # 1 female, 30% community 2 female, 20% community 1 male, 30% community 2 male. table(mydf$sex, mydf$community) # ===================================================================== # SETTING UP THE EMPTY NETWORK # ===================================================================== # The target is a simulated population of 2,000 that carries the sample's sex # and community composition. mynet <- network_initialize(2000) # Sex: 1,000 females (1) and 1,000 males (2). We name the attribute "group" # rather than "sex" because "group" is a special attribute in EpiModel that # enables two-group functionality in the epidemic models (group-specific # infection and recovery rates, group-specific initial conditions, and so on). # See the Module 4 nodal attributes tutorial for the details. sex <- c(rep(1, 1000), rep(2, 1000)) mynet <- set_vertex_attribute(mynet, "group", sex) # Community: within each sex, 400 in community 1 and 600 in community 2. That # reproduces the sample's 40/60 community split within sex, giving 800 in # community 1 and 1,200 in community 2 overall. cmty <- c(rep(1, 400), rep(2, 600), rep(1, 400), rep(2, 600)) mynet <- set_vertex_attribute(mynet, "cmty", cmty) # Confirm the attribute combinations are as expected: 400/600 by 400/600. table(get_vertex_attribute(mynet, "group"), get_vertex_attribute(mynet, "cmty")) # ===================================================================== # THE HINT METRICS, CALCULATED FROM THE DATA # ===================================================================== # The lab suggests several metrics worth calculating before specifying any # terms. The course page states the resulting numbers; here we compute them, # so the derivation is visible rather than taken on trust. Every target # statistic set below traces back to one of these five numbers. # 1. Mean degree by sex. Because every reported tie is cross-sex, each sex's # reports are a complete accounting of that sex's tie-ends. tapply(mydf$num.partners, mydf$sex, mean) # F = 0.8, M = 0.9 # The underlying degree distributions, as proportions within each sex. Nobody # reports more than 2 partners. prop.table(table(mydf$sex, mydf$num.partners), margin = 1) # 2. Mean degree by community. Community 2 has 12 respondents reporting 8 # partnerships in total, a mean degree of 0.667; community 1 has 8 # respondents reporting 9, a mean degree of 1.125. Balance does not bind # here the way it does across sex: ties are not all cross-community (most # of them are within), so there is no cross-community counting identity to # impose. table(mydf$community) tapply(mydf$num.partners, mydf$community, mean) # cmty 1 = 1.125, cmty 2 = 8/12 # 3. Proportion of all ties that are within-community. A tie is within # community when the partner's community matches the respondent's own. We # compare the respondent's community against each partner slot separately # and add the two counts. The NA slots (unreported partners) drop out. n.ties <- sum(mydf$num.partners) # 17 ties reported n.within <- sum(mydf$p1.community == mydf$community, na.rm = TRUE) + sum(mydf$p2.community == mydf$community, na.rm = TRUE) n.ties n.within # 11 within-community n.within / n.ties # 11/17 = 0.647 # 4 and 5. Proportion of each sex reporting exactly one tie. These feed the # sex-specific monogamy terms at the end. tapply(mydf$num.partners == 1, mydf$sex, mean) # F = 0.6, M = 0.3 # ===================================================================== # SPECIFYING THE MODEL, ONE TERM AT A TIME # ===================================================================== # We work iteratively: decide on one term, calculate its target statistic, # and carry forward what that decision has already fixed. Several of these # decisions have no single correct answer. What matters is being explicit # about the assumption each choice implies. # ---- Overall edges ---- # The term is ~edges. To set its target we must reconcile four facts: # # 1. In the observed data, females and males do NOT report the same mean # degree: 0.8 for females, 0.9 for males. # 2. In the model population, the total edges females have with males must # equal the total edges males have with females. This is balance, and it # holds by construction in any actual network. # 3. The two groups are equal in size, both in the sample (10 and 10) and in # the simulated population (1,000 and 1,000). # 4. All ties are cross-sex. # # Points 4 and 2 together say the total tie-ends held by males must equal the # total held by females. Adding point 3, equal group sizes, means that IN THIS # CASE mean degree must also be equal across the sexes. # # That last step depends on the equal group sizes. With unequal groups, # balanced tie counts imply UNEQUAL mean degrees: if 200 females at mean # degree 0.8 hold 160 tie-ends, then 180 males must also hold 160 tie-ends, a # mean degree of 160/180 = 0.889, not 0.8. # # Computing the edge count from each sex's own reports gives two different # answers, which is the problem balance forces us to confront: 1000 * 0.9 # num males * mean degree of males (all ties are F-M) 1000 * 0.8 # num females * mean degree of females (all ties are F-M) # The two are not equal, though they are not far apart, and the gap is easily # within what sampling variation could produce in a sample of 20. Our options: # # 1. Treat the female reports as accurate and the male reports as biased # (why might males over-report?). # 2. Treat the male reports as accurate and the female reports as biased # (why might females under-report?). # 3. Treat both as equally accurate and place the true value in between. # 4. Revise the assumption about the sex ratio itself, so that unequal mean # degrees become consistent with balance. # # We take option 3, splitting the difference at a mean degree of 0.85 for both # sexes. Nothing in the data adjudicates between these options; option 3 simply # declines to privilege either sex's reports. (This is also what ergm.ego does # by default: it treats the discrepancy as sampling error and averages the two # sides' reported tie-ends.) 2000 * 0.85 / 2 # num people * mean degree / 2, so 850 edges # ---- Cross-sex constraint ---- # The population is defined as purely heterosexual, so we need to prohibit # within-sex ties. The term is nodematch("group") with a target statistic of # 0: the count of ties joining two nodes of the same group is held at zero, # which forces every tie to be cross-sex. # # Because 0 is the smallest attainable value of that statistic, its # coefficient is fitted at -Inf, the exact representation of "this never # happens" rather than "this is rare". netest() reports that as it fits, and # the message is expected here rather than a sign of trouble. # # Without this term the model would form female-female and male-male ties, # which would break both the substantive premise of the lab and the balance # logic we just used to set the edges target. # ---- Degree 2 constraint ---- # Nobody in the data reports more than 2 current partners. We carry that into # the model with degrange(from = 3), the count of nodes with degree 3 or more, # again with a target statistic of 0. # # What this assumption costs: a hard zero says degree 3 is IMPOSSIBLE, not # merely unobserved. With 20 respondents, seeing no one at degree 3 is weak # evidence that no one exists at degree 3, and a hard ceiling removes exactly # the high-degree nodes that tend to matter most for transmission. We could # instead assign a small positive target, allowing the occasional higher # degree while keeping it rare. This is another reason to want a much larger # sample before fitting a model of this kind. # ---- Mean degree by community ---- # The two communities differ in size and appear to differ in mean degree, so # we add a nodefactor("cmty") term. # # Why only ONE statistic for a two-category attribute: nodefactor omits a # reference category by default (the value that comes first numerically or # alphabetically, here community 1), and reports the total degree of nodes in # each remaining category. Community 1's total degree is already pinned by the # edges term once community 2's is set, so a second statistic would be # redundant. We therefore need only the expected total degree of community 2. # # From the data: 12 respondents in community 2 report 8 partnerships between # them, a mean degree of 8/12 = 0.667. Applied to the 1,200 community 2 # members of the simulated population: (8/12) * 1200 # mean degree for cmty 2 * pop size of cmty 2, so 800 # Worth a sanity check on what that implies for the other community. There are # 2 * 850 = 1,700 tie-ends in total, so community 1's 800 members hold the # remaining 900, an implied mean degree of 1.125. That matches the sample. # ---- Mixing by community ---- # The data show a preference for within-community ties, so we want an # assortative mixing term. nodematch on a categorical attribute comes in two # flavors: differential (diff = TRUE) adds one term and statistic per value of # the attribute, and uniform (diff = FALSE) adds a single term representing a # global tendency. Which one can we afford? # # THE DEGREES OF FREEDOM ARGUMENT: # # 1. By community, there are only three possible tie types: (1,1), (1,2) and # (2,2). (2,1) is the same as (1,2) because ties are undirected. Counting # the ties of each type starts us with 3 degrees of freedom. # 2. Those three counts must sum to the edges target we already set. That # removes one, leaving 2 d.f. # 3. The number of (1,2) ties plus twice the number of (2,2) ties must equal # the community 2 total degree we set with nodefactor. That removes # another, leaving 1 d.f. # # With one degree of freedom left there is no choice to make: we must use # nodematch("cmty", diff = FALSE). Specifying diff = TRUE would ask for two # statistics where only one is free, over-parameterizing the model. Put # another way, given the overall and community-specific mean degrees already # set, fixing the homophily level for one community automatically fixes it for # the other. # # The statistic is the total number of within-community ties. From the data, # 11 of the 17 reported ties are within-community, so: (11/17) * 850 # proportion within-community * total edges, so 550 # ---- Sex-specific monogamy terms ---- # Finally, the sex-specific tendency toward having exactly one partner, using # degree(1, by = "group"). This yields two target statistics: the expected # number of females and the expected number of males with exactly 1 tie. # # The naive approach reads the proportions straight off the data, where 6/10 # females and 3/10 males report exactly one tie: (6/10) * 1000 # females, 600 (3/10) * 1000 # males, 300 # WE REJECT THOSE TARGETS, for the following reason. We already changed the # sex-specific mean degree relative to the data when we settled the edges # term: the data implied 0.8 for females and 0.9 for males, and balance forced # us to 0.85 for each. If the mean degree moved, the degree distribution that # generates it has to move somewhere too. Holding the degree 1 proportions # fixed at the observed values pushes all of the adjustment onto degrees 0 and # 2, implying these distributions: # # Degree 0 1 2 mean # Females 27.5% 60.0% 12.5% 0.85 # Males 42.5% 30.0% 27.5% 0.85 # # (Observed, for comparison: Females 30.0/60.0/10.0 with mean 0.8; Males # 40.0/30.0/30.0 with mean 0.9.) # # That may be what you want. Or you may want a different assumption. Imagine # you had a larger data set that showed this difference in mean degree, so it # was not so easily attributed to random sampling variation, and suppose you # had some other reason to think that the slight under-reporting on the female # side came from women with 2 partners saying they had 1, and the slight # over-reporting on the male side came from men with 0 partners saying they # had 1. Then you would want to model the degree distributions as: # # Degree 0 1 2 mean # Females 30.0% 55.0% 15.0% 0.85 # Males 45.0% 25.0% 30.0% 0.85 # # We proceed with this second option. It rests on an assumption you may # reasonably dispute. Note, though, that taking the degree 1 statistics # straight from the data is equally an assumption, about how the OTHER degrees # absorbed the change; it is just more hidden. Running the epidemic under both # parameterizations as a sensitivity analysis, and checking whether the # outcomes of interest differ meaningfully, is the cleaner way to handle the # ambiguity. 0.55 * 1000 # females at degree 1, so 550 0.25 * 1000 # males at degree 1, so 250 # ===================================================================== # FINAL FORMATION FORMULA AND TARGET STATISTICS # ===================================================================== # The seven statistics, in formula order: edges (850); same-sex ties (0); # nodes with degree 3+ (0); total degree of community 2 (800); within-community # ties (550); females with degree 1 (550); males with degree 1 (250). formation <- ~edges + nodematch("group", diff = FALSE) + degrange(from = 3) + nodefactor("cmty") + nodematch("cmty", diff = FALSE) + degree(1, by = "group") target.stats <- c(850, 0, 0, 800, 550, 550, 250) # Fit the TERGM. The dissolution model uses a mean relational duration of 187 # time steps. That number does not come from the sample; we never asked the # respondents about duration, so it is taken from the literature on the # closest population that has been surveyed. Not ideal, and a reviewer would # be within their rights to object. Note that it costs no degrees of freedom # on the formation side: it is an offset, a value we impose rather than # estimate. myfit <- netest(mynet, formation = formation, target.stats = target.stats, coef.diss = dissolution_coefs(~offset(edges), 187)) # ===================================================================== # NETWORK DIAGNOSTICS # ===================================================================== # Simulate from the fitted model and compare the realized statistics against # the targets. This checks that the seven constraints we specified are # mutually satisfiable and that the fitting procedure found them. mydx <- netdx(myfit, nsims = 25, nsteps = 500, ncores = 5) mydx # The same information as a plot: simulated statistics against target lines. plot(mydx) # The raw statistics behind the numerical summaries and plots are available # if you want to check them directly or plot them yourself. head(get_nwstats(mydx), 10) # ===================================================================== # EPIDEMIC SIMULATION # ===================================================================== # An SIS epidemic on the fitted network, with the parameters the lab specifies: # a per-act transmission probability of 0.6, an average of 1.8 acts per # partnership per time step, and a recovery rate of 0.1 (a mean infectious # duration of 10 time steps). The .g2 arguments set the corresponding values # for group 2 (males); here they match group 1. myparam <- param.net(inf.prob = 0.6, inf.prob.g2 = 0.6, act.rate = 1.8, rec.rate = 0.1, rec.rate.g2 = 0.1) # Initial conditions: 100 infected females and 100 infected males. myinit <- init.net(i.num = 100, i.num.g2 = 100) # Controls. We pass nwstats.formula to monitor network statistics during the # epidemic simulation, including degrees 0 through 5 by group, to confirm that # no nodes reach a degree above 2 once the epidemic dynamics are running. mycontrol <- control.net("SIS", nsteps = 50, nsims = 5, nwstats.formula = ~edges + nodematch("cmty") + degree(0:5, by = "group"), verbose = TRUE) mySIS <- netsim(myfit, param = myparam, init = myinit, control = mycontrol) # Prevalence by sex. Prevalence in group 1 (females) runs higher than in group # 2 (males). This would once again seem to be driven by the sex-specific # degree distributions we chose. Note that they differ in two ways at once: # males hold more of the concurrency (30% at degree 2 versus 15%), but also # more of the isolates (45% at degree 0 versus 30%), and a node with no # partners cannot be infected. Which of the two dominates here is a conjecture # worth testing rather than a result we have established; re-running with each # difference removed in turn would settle it. plot(mySIS, y = c("i.num", "i.num.g2"), legend = TRUE) # Numerical summaries of sex-specific prevalence at the final time step. mySIS <- mutate_epi(mySIS, prev.f = i.num / num, prev.m = i.num.g2 / num.g2) df <- as.data.frame(mySIS, out = "mean") df$prev.f[df$time == 50] df$prev.m[df$time == 50] # Finally, the network statistics as a time series from within the epidemic # simulation. These should hold on their targets throughout. If they drift, # the epidemic dynamics have had an unintended effect on network structure. plot(mySIS, type = "formation", qnts = FALSE, sim.lines = TRUE)