---
title: "Valued ERGMs, in depth"
engine: knitr
editor: source
subtitle: "Ancillary self-study · ICPSR Network Analysis: Advanced Topics"
---

```{r}
#| label: setup
#| include: false
source("_helpers.R")
require_pkgs(c("ergm", "ergm.count", "network"))
set.seed(6886)
```

{{< include _ancillary_header.qmd >}}

::: {.callout-note title="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$.

```{r}
#| label: data
suppressMessages({
  library(ergm)
  library(ergm.count)
  library(network)
})

data(samplk)
tot <- as.matrix(samplk1) + as.matrix(samplk2) + as.matrix(samplk3)
table(tot)

nw <- network(tot, directed = TRUE, matrix.type = "adjacency",
              ignore.eval = FALSE, names.eval = "nom")
nw
```

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.

```{r}
#| label: checkpoint-data
checkpoint(nodes = network.size(nw),
           nonzero_dyads = sum(tot > 0),
           max_count = max(tot))
```

# 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 |

::: {.callout-warning title="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.
:::

::: {.callout-note collapse="true" title="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.

```{r}
#| label: fits-basic
#| output: false
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)))
```

```{r}
#| label: fits-basic-show
rbind(`sum only`        = c(sum = coef(v1)[["sum"]], nonzero = NA),
      `sum + nonzero`   = c(sum = coef(v2)[["sum"]], nonzero = coef(v2)[["nonzero"]])) |>
  round(4)
```

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.

::: {.callout-important title="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.

```{r}
#| label: fit-mutual
#| output: false
v3 <- fit_or_load("val_mutual",
  ergm(nw ~ sum + nonzero + mutual("min"), response = "nom", reference = ~Poisson,
       control = control.ergm(seed = 6886)))
```

```{r}
#| label: fit-mutual-show
round(coef(v3), 4)
```

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)
```

1. **`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"`.
2. **`combine`**: how to aggregate across all the intermediaries $k$. `"max"` (the single strongest path is what matters) or `"sum"` (they accumulate).
3. **`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.

```{r}
#| label: fit-trans
#| output: false
v4 <- fit_or_load("val_trans",
  ergm(nw ~ sum + nonzero + transitiveweights("min", "max", "min"),
       response = "nom", reference = ~Poisson,
       control = control.ergm(seed = 6886)))
```

```{r}
#| label: fit-trans-show
round(coef(v4), 4)
```

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

::: {.callout-warning title="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.

```{r}
#| label: fit-covar
#| output: false
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)))
```

```{r}
#| label: fit-covar-show
round(coef(v5), 4)
```

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.

::: {.callout-note collapse="true" title="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:

1. **`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.
2. **`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).

```{r}
#| label: mcmc-diag
#| fig-width: 8
#| fig-height: 5
#| results: hide
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.

```{r}
#| label: sim
#| fig-width: 8
#| fig-height: 3.4
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)
}
```

```{r}
#| label: sim-table
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)
```

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.

::: {.callout-tip collapse="true" title="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.

```{r}
#| label: sim-oom
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))
```

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

::: {.panel-tabset}

## Core

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.

```{r}
#| eval: false
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?

## Stretch

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.

## Solution

```{r}
#| label: ex-binom
#| output: false
v_binom <- fit_or_load("val_binom",
  ergm(nw ~ sum + nonzero, response = "nom", reference = ~Binomial(3),
       control = control.ergm(seed = 6886)))
```

```{r}
#| label: ex-binom-show
rbind(Poisson  = round(coef(v2), 4),
      Binomial = round(coef(v_binom), 4))
```

**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.

```{r}
#| label: ex-trans2
#| output: false
v4b <- fit_or_load("val_trans_geomean",
  ergm(nw ~ sum + nonzero + transitiveweights("geomean", "sum", "geomean"),
       response = "nom", reference = ~Poisson,
       control = control.ergm(seed = 6886)))
```

```{r}
#| label: ex-trans2-show
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)
AIC(v4, v4b)
```

**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:

1. **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?
2. **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.
3. **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.

```{r}
#| label: versions
#| cache: false
#| code-fold: true
#| code-summary: "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)
cat("R", as.character(getRversion()), "\n")
```
