Valued ERGMs, in depth

Ancillary self-study · ICPSR Network Analysis: Advanced Topics

Author

Shahryar Minhas

Published

July 22, 2026

Optional: not covered in the 2026 course

This document is not part of any taught session. Nobody is going to lecture it at you, nobody is going to grade you on it, and you can ignore it completely without falling behind.

It exists because it is a thing you will plausibly need next, and because it is easier to hand you something that works than to hand you a citation and wish you luck. Everything here runs end to end on the same machine setup you used in class. Every number and every figure on this page was produced by a chunk in this file: there are no screenshots of R output anywhere in the ancillary package.

How to work through it alone. Get the source (click the </> Code button at the top right, then View Source, and copy it into a .qmd file: or open the .qmd you were given directly), run the chunks in order, and read the prose between them: the prose is where the argument lives. Slow fits are identified in the prose, and precomputed fits load a shipped result while still showing you the code that produced it. Folded Depth boxes are optional extras; folded Stuck? boxes are hints. Open whichever you need and ignore the other. Exercises have a Solution tab: you are working alone, so use it.

What this extends, and why it is here

Extends: the ERGM session (day 11, W3 Mon), block C3. There we introduced valued ERGMs: why you need a reference measure, sum as the valued analogue of edges, and the sum + nonzero sign flip. That is the whole idea and it is genuinely enough to get started.

This document is the rest of it. Specifically: the other reference measures and how to choose between them, valued mutuality, the transitiveweights family (valued triadic closure, which is much stranger than it sounds), valued covariate effects, and how to simulate from a valued fit and check it. If you have valued data: counts of interactions, message volumes, trade values, durations: this is where the terms you actually need live.

What you need: ergm, ergm.count, network. All on CRAN, all installed for the course.

1. Where we left off

The monastery again. Sampson observed his novices at three time points, so summing the three adjacency matrices gives a count from 0 to 3 for each ordered pair: how many times \(i\) nominated \(j\).

suppressMessages({
  library(ergm)
  library(ergm.count)
  library(network)
})

data(samplk)
tot <- as.matrix(samplk1) + as.matrix(samplk2) + as.matrix(samplk3)
table(tot)
#> tot
#>   0   1   2   3 
#> 236  38  20  30

nw <- network(tot, directed = TRUE, matrix.type = "adjacency",
              ignore.eval = FALSE, names.eval = "nom")
nw
#>  Network attributes:
#>   vertices = 18 
#>   directed = TRUE 
#>   hyper = FALSE 
#>   loops = FALSE 
#>   multiple = FALSE 
#>   bipartite = FALSE 
#>   total edges= 88 
#>     missing edges= 0 
#>     non-missing edges= 88 
#> 
#>  Vertex attribute names: 
#>     vertex.names 
#> 
#>  Edge attribute names: 
#>     nom

Note the three arguments that make this a valued network rather than a binary one: ignore.eval = FALSE tells network() not to throw the weights away, and names.eval = "nom" gives the edge attribute a name you will pass as response = "nom" to every model below. Leave either out and you will silently fit a binary ERGM to dichotomised data, which is a mistake that does not announce itself.

checkpoint(nodes = network.size(nw),
           nonzero_dyads = sum(tot > 0),
           max_count = max(tot))
#> ------------------------------------------------------------------
#> CHECKPOINT: nodes = 18   |   nonzero_dyads = 88   |   max_count = 3
#> ------------------------------------------------------------------

2. Reference measures: the thing that has no binary analogue

This is the concept the ERGM session had time to name but not to develop, and it is the one that actually distinguishes valued ERGMs from binary ones.

For a binary ERGM, writing down the sufficient statistics is enough. The sample space is fixed: all graphs on \(n\) nodes, and every graph in it starts equally likely. Specify \(g(y)\) and you have specified the model.

For a valued ERGM that is not true, and the reason is worth internalising. Suppose your statistic is sum, the total of all edge values. That pins down the mean. It says nothing about whether the counts should be Poisson-shaped, geometric-shaped, or bounded at some maximum. Infinitely many distributions have the same total. So you must supply a reference measure: the baseline distribution the model tilts away from, which fixes both the sample space and the shape.

\[\Pr(Y = y) = \frac{h(y)\exp\{\theta^{\top}g(y)\}}{\kappa(\theta)}\]

where \(h(y)\) is the reference. Setting \(h(y) = 1\) over binary graphs recovers the ordinary ERGM, so the binary case is the special case where you did not have to think about this.

reference	\(h(y)\)	sample space	use it when
`~Poisson`	\(\prod 1/y_{ij}!\)	counts \(0, 1, 2, \ldots\)	counts with no natural ceiling: messages, citations, interactions
`~Binomial(trials)`	\(\prod \binom{n}{y_{ij}}\)	\(0, \ldots, n\)	counts out of a known maximum: nominations from a fixed roster
`~Geometric`	\(1\)	counts \(0, 1, 2, \ldots\)	heavier tails than Poisson, but read the warning below
`~DiscUnif(a,b)`	\(1\) on \([a,b]\)	integers \(a\) to \(b\)	bounded ratings: a 1–5 Likert tie strength

The geometric reference is a degeneracy trap

~Geometric has \(h(y) = 1\), which means it puts no penalty at all on large values. Combine that with a positive sum coefficient and the model happily wanders off toward infinite edge weights. The ergm.count vignette is explicit that the geometric reference is prone to degeneracy for exactly this reason and that Poisson is the safer default for counts.

If you think you need heavier tails than Poisson, the better move is usually to keep the Poisson reference and add a dispersion-handling term, rather than to change the reference.

Depth: why the reference is not just a prior

It is tempting to file the reference measure under “prior” and move on. Resist that, because it does something a prior does not: it defines the sample space. Under ~Poisson the space is all non-negative integer-valued arrays; under ~DiscUnif(1,5) it is arrays with entries in \(\{1,\ldots,5\}\) and there is no such thing as a zero. Those are different models of the world, not different beliefs about the same one.

The practical consequence you will hit immediately: likelihood-based comparisons are only valid within a reference. ergm prints this warning on every valued fit:

“Null model likelihood calculation is not implemented for valued ERGMs at this time. This means that all likelihood-based inference (LRT, Analysis of Deviance, AIC, BIC, etc.) is only valid between models with the same reference distribution and constraints.”

So you cannot use AIC to choose between a Poisson model and a geometric one. You choose the reference from what you know about how the data were generated, and then you use AIC within it. This is a case where the software’s limitation and the correct statistical practice happen to coincide.

3. The `sum + nonzero` result, restated

The taught session’s headline, reproduced here so this document stands alone.

v1 <- fit_or_load("val_sum",
  ergm(nw ~ sum, response = "nom", reference = ~Poisson,
       control = control.ergm(seed = 6886)))

v2 <- fit_or_load("val_sum_nonzero",
  ergm(nw ~ sum + nonzero, response = "nom", reference = ~Poisson,
       control = control.ergm(seed = 6886)))

rbind(`sum only`        = c(sum = coef(v1)[["sum"]], nonzero = NA),
      `sum + nonzero`   = c(sum = coef(v2)[["sum"]], nonzero = coef(v2)[["nonzero"]])) |>
  round(4)
#>                   sum nonzero
#> sum only      -0.6023      NA
#> sum + nonzero  0.3778 -2.1068

The sum coefficient flips sign. On its own it is negative; add nonzero and it turns positive, while nonzero itself comes out strongly negative.

The reading: a bare Poisson reference cannot accommodate how many zeros are in real network data. Left to explain 236 zero cells with a single parameter, sum is dragged negative: it is doing the job of a zero-inflation term badly. Give the model a dedicated zero-inflation parameter and sum is freed to describe the counts that actually exist. Conditional on a tie existing at all, the monks nominate each other more than Poisson predicts, not less. The original negative was zero-inflation wearing a disguise.

The general lesson, which is not about monks

Almost all valued network data is zero-inflated, because most pairs of anything never interact. nonzero (or Binary(), or a hurdle construction) is close to mandatory, not optional.

If you fit a valued ERGM with only sum and report the coefficient, you are very likely reporting a zero-inflation parameter and calling it a volume parameter. Check by adding nonzero and seeing whether anything moves. If it moves a lot, you had the wrong model.

4. Beyond `sum`: the terms you will actually need

4.1 Mutuality, valued

Binary mutual asks whether \(i \to j\) and \(j \to i\) co-occur. Valued mutuality has to answer a harder question: what does it mean for a weight of 3 and a weight of 1 to be “reciprocated”?

ergm.count gives you a choice of how much correlation to credit, via the form argument:

mutual("min"): credits \(\min(y_{ij}, y_{ji})\). Conservative: a pair is reciprocal only up to the weaker of the two directions.
mutual("nabsdiff"): credits \(-|y_{ij} - y_{ji}|\). Penalises imbalance rather than rewarding overlap.
mutual("product"): credits \(y_{ij} \cdot y_{ji}\), i.e. raw correlation.

v3 <- fit_or_load("val_mutual",
  ergm(nw ~ sum + nonzero + mutual("min"), response = "nom", reference = ~Poisson,
       control = control.ergm(seed = 6886)))

round(coef(v3), 4)
#>        sum    nonzero mutual.min 
#>     0.0164    -2.2077     1.4210

Strong positive mutuality, and notice what it does to sum: adding mutuality pulls sum back toward zero, because some of the volume the model was attributing to a general propensity to nominate is in fact a propensity to nominate back.

4.2 Triadic closure, valued, and why it is genuinely hard

Here is the part of valued ERGMs that catches people out.

In a binary network, “\(i\) and \(j\) share a partner \(k\)” is a yes/no question. In a valued network it is not. If \(y_{ik} = 3\) and \(y_{kj} = 1\), how strong is that two-path? And if \(i\) has ten weak two-paths to \(j\), does that count for more or less than one strong one?

There is no single right answer, so ergm.count makes you choose three things, which is why transitiveweights() takes three arguments:

transitiveweights(twopath, combine, affect)

twopath: how to combine \(y_{ik}\) and \(y_{kj}\) into the strength of one two-path. "min" (a chain is as strong as its weakest link) or "geomean".
combine: how to aggregate across all the intermediaries \(k\). "max" (the single strongest path is what matters) or "sum" (they accumulate).
affect: how the aggregated two-path strength relates to \(y_{ij}\) itself. "min" or "geomean".

The default, transitiveweights("min", "max", "min"), reads as: a two-path is as strong as its weakest leg; what matters is the strongest available two-path; and closure is credited up to the weaker of that path and the direct tie. That is a defensible reading of “triadic closure” for valued data. It is not the only one.

v4 <- fit_or_load("val_trans",
  ergm(nw ~ sum + nonzero + transitiveweights("min", "max", "min"),
       response = "nom", reference = ~Poisson,
       control = control.ergm(seed = 6886)))

round(coef(v4), 4)
#>                           sum                       nonzero 
#>                        0.2455                       -2.3995 
#> transitiveweights.min.max.min 
#>                        0.2656

Positive, as expected: monks nominate those their nominees nominate.

Report your three arguments. Always.

Because transitiveweights has three switches, transitiveweights is not a term, it is a family of terms. Two papers both reporting “a positive valued transitivity effect” may have fitted quite different models.

So: state all three arguments in your write-up, and if you had a substantive reason for choosing them, give it. If you did not: if you took the default because it was the default: say that instead of inventing a rationale after the fact. The defaults are sensible and there is nothing wrong with using them.

There is also a practical warning to expect. On the fits above ergm reports Best valid proposal 'DiscTNT' cannot take into account hint(s) 'triadic': the triadic MCMC hint that speeds up binary models with closure terms has no valued counterpart, so it is ignored. This is informational, not an error, but it does mean valued triadic models mix more slowly than their binary equivalents. Budget more iterations and check mcmc.diagnostics().

4.3 Covariate effects: how a nodal attribute maps to volume

Sampson recorded which faction each novice belonged to: Loyal Opposition, Young Turks, Outcasts, Waverers. A natural valued question: do monks in the same faction nominate each other more heavily, not merely more often?

That last distinction is the point. Binary nodematch("group") can only ask whether same-faction ties are more likely to exist. The valued version, nodematch("group", form = "sum"), asks whether same-faction ties carry more weight: a question that only exists once edges have values.

network::set.vertex.attribute(nw, "group", as.character(samplk1 %v% "group"))

v5 <- fit_or_load("val_nodematch",
  ergm(nw ~ sum + nonzero + nodematch("group", form = "sum"),
       response = "nom", reference = ~Poisson,
       control = control.ergm(seed = 6886)))

round(coef(v5), 4)
#>                 sum             nonzero nodematch.sum.group 
#>             -0.2924             -1.5687              1.0275

Strongly positive: conditional on the reference and the zero-inflation term, monks nominate same-faction brothers with noticeably higher counts. The form = "sum" argument is doing the work: it tells ergm to accumulate the weights on matching dyads rather than count matching edges. Drop it and you get the binary-style match term applied to a valued network, which is rarely what you want and does not warn you.

Depth: the dispersion terms, and why they are not in this document

There is a second family of valued nodal terms: nodeocovar, nodeicovar, nodecovar: that measure dispersion rather than level: not “do high-attr actors send more” but “do actors differ in how evenly they spread what they send.” These are conceptually the valued analogue of degree terms, and on the right data they are exactly what you want.

They are deliberately absent from the runnable part of this document, for two honest reasons I hit while building it:

nodecovar requires an undirected network. On directed Sampson it errors outright (Term may not be used with networks with directed==TRUE), and nodesqrtcovar is deprecated in favour of nodecovar(transform = "sqrt"), which inherits the same restriction.
nodeocovar / nodeicovar do run on directed data but mix very slowly here: repeated fits produced no reliable convergence on 18 nodes, because these are strongly dyad-dependent terms and the triadic MCMC hint is unavailable for valued models (§4.2).

The lesson is more useful than a fitted table would have been: valued dyad-dependent terms are substantially harder to estimate than their binary cousins, and not every term that exists is practical on every dataset. If you need dispersion terms, use an undirected network, budget real MCMC time, and read mcmc.diagnostics() carefully before trusting anything. The ergm.count vignette works them on undirected data for exactly this reason.

5. Checking a valued fit

Two things to do, and neither is optional.

5.1 MCMC diagnostics

Every model in §4 is dyad-dependent and was fitted by MCMC, so the chains need looking at. The valued case deserves more attention than the binary case, not less, because the triadic proposal hint is unavailable (§4.2).

mcmc.diagnostics(v3, which = "plots")

What you want: traces that look like fuzzy horizontal bands with no drift, and densities that look roughly unimodal. What you do not want: a trace that wanders steadily in one direction, or one that sticks at a value for long stretches.

5.2 Simulate and compare

The valued analogue of gof(). As of ergm 4.11, gof() does support valued models, but the transparent thing, and the thing that generalises to any statistic you care about: is to simulate and compare by hand.

set.seed(6886)
sims <- simulate(v3, nsim = 200, response = "nom", output = "network")

obs_stats <- summary(nw ~ sum + nonzero + mutual("min"), response = "nom")
sim_stats <- t(sapply(sims, function(s)
  summary(s ~ sum + nonzero + mutual("min"), response = "nom")))

par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
for (j in seq_along(obs_stats)) {
  hist(sim_stats[, j], breaks = 20, col = "grey85", border = "white",
       main = names(obs_stats)[j], xlab = "simulated value")
  abline(v = obs_stats[j], col = "firebrick", lwd = 3)
}

data.frame(
  statistic = names(obs_stats),
  observed  = round(obs_stats, 2),
  sim_mean  = round(colMeans(sim_stats), 2),
  p_value   = round(sapply(seq_along(obs_stats), function(j) {
    p <- mean(sim_stats[, j] >= obs_stats[j])
    2 * min(p, 1 - p)
  }), 3),
  row.names = NULL)
#>    statistic observed sim_mean p_value
#> 1        sum      168   168.56    0.97
#> 2    nonzero       88    87.60    0.98
#> 3 mutual.min       47    47.63    0.98

The red line is the observed value; the histogram is what the model produces. For in-model statistics this should land comfortably inside the distribution: the model is fitted to reproduce exactly these, so if it does not, something has gone wrong with estimation rather than with the model.

Stuck? abline(v = obs) drew three lines and I wanted one

A classic. obs_stats is a vector of three numbers, and abline() is vectorised: pass it a vector and it draws one line per element, on whichever panel is current.

Inside the loop above the fix is abline(v = obs_stats[j]), indexing to the one you want. If you find yourself with a plot carrying mysterious extra vertical lines, this is almost always why.

The real test is an out-of-model statistic: something the model was not fitted to reproduce. Model v3 has sum + nonzero + mutual("min") in it, so valued transitivity is fair game: the model was never told about triadic closure, so whether it reproduces the observed amount is a genuine test.

set.seed(6886)
obs_oom <- summary(nw ~ transitiveweights("min", "max", "min"), response = "nom")
sim_oom <- sapply(sims, function(s)
  summary(s ~ transitiveweights("min", "max", "min"), response = "nom"))

c(observed = round(obs_oom, 2),
  sim_mean = round(mean(sim_oom), 2),
  sim_q025 = round(quantile(sim_oom, 0.025), 2),
  sim_q975 = round(quantile(sim_oom, 0.975), 2))
#> observed.transitiveweights.min.max.min                               sim_mean 
#>                                  96.00                                  85.40 
#>                          sim_q025.2.5%                         sim_q975.97.5% 
#>                                  47.98                                 139.13

Compare the observed value against the simulated interval. Here the observed transitivity sits comfortably inside the simulated band: the model reproduces a feature it was never told about, which is genuine (if modest) evidence in its favour: the mutuality and volume terms it does contain already generate about the right amount of closure, so you would not gain much by adding a transitiveweights term. Had the observed value landed above the band instead, you would have found closure the model is missing, and §4.2 has the term to add. That is the whole logic of an out-of-model check: the informative outcome is the one where observed and simulated disagree.

6. Exercises

Fit the sum + nonzero model with a binomial reference instead of Poisson. The counts run 0 to 3 because there were three observation waves, so trials = 3 is the natural choice.

v_binom <- ergm(nw ~ sum + nonzero,
                response  = "nom",
                reference = ~Binomial(______),
                control   = control.ergm(seed = 6886))

Then answer: the coefficients differ from the Poisson version. In one sentence, why can you not use AIC to decide which of the two is better?

Refit §4.2’s transitivity model with transitiveweights("geomean", "sum", "geomean") instead of the default ("min", "max", "min").

Report both coefficients. Then make the harder argument: which of the two specifications is a better description of this dataset, and what feature of the data would you look at to decide? “The one with the bigger coefficient” is not an answer. You may compare AIC because the two models share a response, reference, and constraints, but explain why AIC alone cannot tell you which closure interpretation is substantively right.

v_binom <- fit_or_load("val_binom",
  ergm(nw ~ sum + nonzero, response = "nom", reference = ~Binomial(3),
       control = control.ergm(seed = 6886)))

rbind(Poisson  = round(coef(v2), 4),
      Binomial = round(coef(v_binom), 4))
#>             sum nonzero
#> Poisson  0.3778 -2.1068
#> Binomial 0.3956 -3.5603

Core. The sum coefficients are close (both positive, around 0.38–0.40), but look at nonzero: it is markedly more negative under the binomial reference (around −3.6) than under Poisson (around −2.1). That is the reference doing visible work. A binomial with trials = 3 already knows the counts are capped at 3 and puts more of its own mass on the low end, so relative to that stricter baseline the excess of zeros looks even more extreme, and nonzero has to push harder. The sign-flip story from §3 survives the change of reference, but the magnitudes do not transfer, which is the whole point of the next paragraph.

You cannot use AIC to choose between them because they are not defined over the same sample space. The Poisson model assigns probability to a monastery where someone is nominated 47 times; the binomial model assigns that outcome probability zero, because it does not exist in its universe. Likelihoods computed over different sample spaces are not comparable numbers, and ergm warns you about exactly this on every valued fit (see the Depth box in §2).

Choose the reference from what you know about the data-generating process. Here the binomial is arguably the more honest choice, since three waves genuinely cannot produce a count above three.

v4b <- fit_or_load("val_trans_geomean",
  ergm(nw ~ sum + nonzero + transitiveweights("geomean", "sum", "geomean"),
       response = "nom", reference = ~Poisson,
       control = control.ergm(seed = 6886)))

data.frame(
  spec = c("min / max / min", "geomean / sum / geomean"),
  transitivity = c(round(coef(v4)[3], 4), round(coef(v4b)[3], 4)),
  sum          = c(round(coef(v4)[1], 4), round(coef(v4b)[1], 4)),
  row.names    = NULL)
#>                      spec transitivity    sum
#> 1         min / max / min       0.2656 0.2455
#> 2 geomean / sum / geomean       0.2040 0.1806
AIC(v4, v4b)
#>     df       AIC
#> v4   3 -160.1729
#> v4b  3 -162.5537

Stretch. The two coefficients are on different scales: they are not measuring the same quantity in different amounts, they are measuring different quantities. ("min","max","min") counts the single strongest two-path and caps its contribution at the weaker leg; ("geomean","sum","geomean") accumulates across all intermediaries and lets many weak paths add up. Naturally the summing version produces a larger raw number. Comparing the magnitudes is meaningless.

AIC is legitimate here. The models share a response, sample space, reference, and constraints, and AIC does not require models to be nested. A lower value says one specification has the better expected information tradeoff under the fitted likelihood. It does not tell you that its definition of a valued two-path is the true social mechanism. Choosing the statistic still requires theory and diagnostics, rather than treating an information criterion as a mechanism identifier.

What you would actually look at, in order:

The distribution of the counts. With a maximum of 3 and most non-zero cells at 1, the difference between “strongest path” and “sum of paths” is small, because there are not many multi-path pairs with much weight to accumulate. On denser valued data, such as email volumes or trade, the choice matters much more. Diagnostic: how many pairs have two or more non-trivial two-paths between them?
What you think closure means here. Sampson’s monks are eighteen people in one building. The plausible mechanism is “I nominate the person my close friend nominates,” which is a strongest-path story, not an accumulation story. That argues for the default.
Out-of-model fit, as in §5.2. Simulate from each and compare on a statistic neither was fitted to, such as the valued triad census or nodeicovar. This complements AIC by asking whether the model reproduces features it was not tuned on.

Item 2 is the one that should carry the most weight, and it is the one people skip.

7. Readings

Krivitsky (2012), “Exponential-family random graph models for valued networks,” Electronic Journal of Statistics 6:1100–1128. The paper the ERGM session assigns. Sections 2 and 3 are the reference-measure argument; section 4 is the term taxonomy in §4 above.
The ergm.count vignette: vignette("valued", package = "ergm.count"). Runnable, current, and maintained by Krivitsky himself. If you read one thing, read this.
Krivitsky, Hunter, Morris & Klumb (2023), “ergm 4: New Features,” Journal of Statistical Software 105(6). Documents valued gof() support and the current term operators.

Versions: expand if your numbers differ from mine

pk <- c("ergm", "ergm.count", "network")
data.frame(package = pk,
           version = sapply(pk, function(p) as.character(packageVersion(p))),
           row.names = NULL)
#>      package version
#> 1       ergm  4.12.0
#> 2 ergm.count   4.1.3
#> 3    network  1.20.0
cat("R", as.character(getRversion()), "\n")
#> R 4.3.3