ergm.triadmatch

An ergm extension providing two terms — gwespmatch and triadmatch — that split edgewise shared partners by the attribute pattern of the triad.

It lets you ask questions that plain ERGM terms cannot: is triadic closure stronger within groups than between them? — and, once a tie does cross a group boundary, whose group is the third party in?

---

The problem it solves

Suppose you find that same-race people are disproportionately tied to each other. The obvious explanation is homophily — a preference for similar others. But there is a second, purely structural explanation: triadic closure. Friends-of-friends become friends, and if your friends are already mostly same-race, closure alone will manufacture same-race ties without anyone preferring anything.

So homophily and closure are confounded. The standard fix is to put both nodematch() and gwesp() in the model and let them compete. That's necessary but not sufficient, because it assumes closure works at the same rate everywhere. If closure is actually stronger inside groups than across them, that interaction is a real mechanism — and leaving it out biases both the homophily and the closure estimates.

gwespmatch is the missing interaction term. Grund & Densley (2015) call it GWESP_MATCH; Hong et al. (2024) call it categorical closure and show that adding it improves model fit and changes the homophily coefficients.

You cannot get this by multiplying nodematch by gwesp. As Grund & Densley put it: "One cannot simply multiply the scores for GWESP and homophily. Instead, one must incorporate a new configuration in the model."

And "within-group closure" is not one configuration but five. A closed triad can be homogeneous, or a within-group tie brokered by an outsider, or a cross-group tie brokered from either side of the boundary. Agneessens, Lomi & Pistocchi (2026) show these carry competing theoretical predictions — contact theory expects the broker on the out-group side, indirect-contact theory on the in-group side — so pooling them can cancel real effects. gwespmatch estimates each separately.

What it counts

For each edge participating in a triad of the chosen type, count that edge's qualifying shared partners, then weight additional partners with geometrically diminishing returns, exactly as gwesp does. This parameterization is designed to reduce the degeneracy risk associated with unweighted triangle counts; it is not a guarantee against degeneracy in every model.

The triad types

A triad has three dyads: the focal tie t–h, the tail-side leg t–k, and the head-side leg k–h. The type is entirely determined by which of those are within-group. And because attribute matching is transitive, you can never have exactly two matching — so only five of the eight patterns exist, the five set partitions of {t, h, k}:

    in.in            in.out           out.tail         out.head         out.none

      ● k              ○ k              ● k              ○ k              ◍ k
     / \              / \              / \              / \              / \
    ●───●            ●───●            ●───○            ●───○            ●───○
    t   h            t   h            t   h            t   h            t   h

  all three       within-group      cross-group      cross-group      cross-group
   match         tie, out-group     tie, third        tie, third      tie, third in
                  third party       on tail's side   on head's side    a third group

● ○ ◍ = three different attribute values.

Read the name as focal.partner: the prefix says whether the focal tie is within-group (in) or cross-group (out), the suffix says whose group the third party shares.

homophily	Condition	Slides
"in.in" (default)	xt = xh = xk	Type A, homogeneous
"in.out"	xt = xh ≠ xk	Type C, close-pair
"out.tail"	xt = xk ≠ xh	Type B, jumping
"out.head"	xt ≠ xk = xh	Type D, alter-jump
"out.none"	all three differ	— (needs ≥ 3 categories)
"dyad"	focal tie matches, k free	pools in.in + in.out
"out.one"	cross-group tie, k matches exactly one endpoint	pools out.tail + out.head

"triad" remains an accepted alias for "in.in", which is still the default and still reproduces the published Grund & Densley / Hong et al. statistic.

out.tail vs out.head is a pure asymmetry of direction — when a tie crosses a group boundary, is the broker on your side or on theirs? On the Sampson monastery data, with decay = 0.5 and the default type = "OTP", those statistics are 13.18 and 8.79. That distinction is not meaningful for undirected networks: ergm stores each undirected edge with tail < head, so either term would silently classify partners by the lower-indexed endpoint. Both are therefore rejected with an error; use their pooled counterpart, "out.one", instead.

Additivity — read this before fitting several types

The five types partition each edge's shared partners, so the counts add up. The GW statistics do not, because 1 − q^ℓ, where q = 1 − e^−α, is applied per edge and is concave in ℓ. Fitting all five is therefore not a decomposition of gwesp, and any reasoning that is linear in triad counts does not survive the GW transform.

If you need the additive version, that is exactly what triadmatch is for.

triadmatch — the unweighted count

triadmatch(match, type = "OTP", homophily = "in.in")

Same triads, same homophily and type vocabulary, no geometric weighting:

g(y) = Σ_edges ESP(edge) = Σ_ℓ ℓ · D_ℓ

i.e. every qualifying shared-partner incidence is counted once for each focal edge. Thus, an undirected triangle contributes 3 (one contribution for each of its three focal edges), whereas a directed transitive triple contributes through the focal edge selected by the chosen shared-partner type. This is the ERGM analogue of the RSiena transitive-triplet effects in the Agneessens et al. slides, which are likewise unweighted — so it is the term to use for a like-for-like comparison with those estimates.

The two terms answer different questions about the same triads:

	gwespmatch	triadmatch
contribution of an edge with ℓ partners	eα(1 − qℓ), concave in ℓ	ℓ
types add up?	no	yes
geometric downweighting	yes — designed to reduce degeneracy risk	no
an edge with 3 partners at decay = 0	scores 1	scores 3

data(sampson)
# in.in + in.out == dyad exactly; the analogous GW statistics need not add.
s <- summary(samplike ~ triadmatch(match = "group", homophily = "in.in")
                        + triadmatch(match = "group", homophily = "in.out")
                        + triadmatch(match = "group", homophily = "dyad"))
stopifnot(unname(s[1] + s[2]) == unname(s[3]))

Reach for triadmatch to compare against triplet counts from other software, or when you want the decomposition to be exact. Reach for gwespmatch when you are fitting a model — unweighted triad counts are what GWESP was invented to avoid.

Installation

# install.packages("remotes")
remotes::install_github("benrosche/ergm.triadmatch")

Requires ergm (>= 4.0) and network.

Getting help

?gwespmatch                 # the weighted term
?triadmatch                 # the unweighted count
?ergm.triadmatch            # package overview and the triad-type table

search.ergmTerms("triad")   # both terms appear in ergm's own term search

Usage

gwespmatch(decay, match, fixed = TRUE, cutoff = 30, type = "OTP",
           homophily = "in.in")

Argument	Default	Description
decay	— (required)	Decay α ≥ 0. Higher α = more credit for additional shared partners. decay = 0 counts each qualifying edge that closes at least one triangle. Fixed if fixed = TRUE; a starting value if fixed = FALSE.
match	— (required)	Name of the vertex attribute defining the groups.
homophily	"in.in"	Triad type — see the table above. "triad" is an alias for the default.
fixed	TRUE	TRUE: one statistic at the given decay. FALSE: curved ERGM, decay estimated jointly with the coefficient.
cutoff	30	Number of ESP bins when fixed = FALSE.
type	"OTP"	Directed shared-partner shape (below). Ignored for undirected networks.

Worked example

library(ergm)
library(ergm.triadmatch)

data(faux.magnolia.high)   # undirected friendship network, with a Race attribute

fit <- ergm(faux.magnolia.high ~ edges
                               + nodematch("Race")                    # dyadic homophily
                               + gwesp(0.5, fixed = TRUE)             # closure in general
                               + gwespmatch(0.5, match = "Race"))     # typed closure within race
summary(fit)

Read the coefficients conditionally, holding the other model statistics fixed:

• nodematch("Race") — association between same-race ties and tie formation at the dyad level.
• gwesp — geometrically weighted closure regardless of race.
• gwespmatch — the additional conditional log-odds contribution associated with its change statistic for closure among three people of the same race.

A positive gwespmatch coefficient is evidence consistent with stronger categorical closure, conditional on the rest of the specification. Compare the nodematch coefficient across models to see how adding typed closure changes the estimated dyadic association, but do not interpret that coefficient change as a causal or additive “share” of homophily.

Directed networks

Choose which shared-partner shape counts, following ergm's convention. For an edge t → h, node k is a shared partner under:

type	Condition on k	Name
"OTP" (default)	t → k → h	outgoing two-path (transitivity)
"ITP"	h → k → t	incoming two-path (cyclicality)
"RTP"	t ↔ k ↔ h	reciprocated two-path
"OSP"	t → k ← h	outgoing shared partner
"ISP"	t ← k → h	incoming shared partner

data(sampson)   # directed liking network among monks, with a "group" attribute
summary(samplike ~ gwespmatch(0.5, type = "OTP", match = "group"))

Estimating the decay (curved ERGM)

fit <- ergm(faux.magnolia.high ~ edges + nodematch("Race")
                               + gwespmatch(0.5, fixed = FALSE, match = "Race"))

decay is then estimated jointly with the coefficient rather than fixed by you.

The statistic

Let q = 1 − e^−α, and let D_ℓ be the number of qualifying edges with exactly ℓ qualifying shared partners — "qualifying" meaning admitted by the chosen type and homophily rule, which for the out.* types means cross-group edges. Then

g(y; α) = e^α · Σ_ℓ (1 − q^ℓ) · D_ℓ

An edge with one shared partner contributes 1 regardless of α; further partners add geometrically less. With fixed = FALSE the term becomes a curved ERGM carrying D₁…D_cutoff and estimating α via the standard GWDECAY map θ_ℓ = θ·e^α(1 − q^ℓ).

See ?gwespmatch for the change statistic and the curved-form gradient.

Validation

Every statistic is checked against:

• Grund & Densley (2015), Table 1 — reproduces their published GWESP_Match values (5.00 / 5.22 / 5.39) exactly.
• ergm itself — gwespmatch(d, match = "a") equals F(~gwesp(d, fixed = TRUE), ~nodematch("a")) on random networks, for all five directed types and undirected.
• An independent pure-R reference — a naive triple loop transcribed straight from the definitions, compared across every triad type × all five directed types × undirected × several decays, on networks with three attribute categories so out.none is reachable.
• The partition identity — the per-edge counts of the five disjoint types sum to the unrestricted shared-partner count, edge by edge.
• The slide figures — one triad at a time, confirming that each configuration fires exactly one type and zero for the other four.
• Change statistics — verified toggle-by-toggle against recomputed summary statistics via ergm.godfather, for every triad type, in the fixed, curved and triadmatch forms. This is also what validates the caches: the change statistics read them, while the summary statistics recount from scratch, so agreement between the two exercises both paths against each other.
• Additivity — triadmatch types are asserted to sum exactly; gwespmatch types are asserted not to, so the documented caveat cannot silently rot.

Run them with:

# From a checkout of this repository:
devtools::test()

Do you even need this package?

For six of the seven gwespmatch homophily settings, yes — there is no equivalent built-in construction.

homophily = "in.in" — often not. ergm's operator syntax already gets you there, and this is exactly how Hong et al. (2024) did it:

gwespmatch(0.5, match = "Race", homophily = "in.in")   # these two
F(~gwesp(0.5, fixed = TRUE), ~nodematch("Race"))       # are the same number

F() restricts the network to within-group edges and then runs gwesp on that subgraph. Since every edge of a homogeneous triad is a within-group edge, the subgraph keeps exactly the triads you want. If you just need the Grund & Densley statistic once, F() avoids installing an additional extension package.

Every other type — no F() equivalent exists. F() filters the network, which is all-or-nothing: it strips a cross-group edge from the two-path legs (t–k and k–h) at the same time as from the focal position. But every other triad type imposes a mixed condition across the three dyads:

• "dyad" and "in.out" need the focal edge filtered while the legs are still evaluated on the full graph — an out-group k must stay visible.
• "out.tail" needs the tail-side leg within-group while the focal tie and the head-side leg are cross-group. There is no single edge predicate that keeps one leg and drops the other, because the same edge plays both roles elsewhere in the network.

No composition of ergm's operators expresses these. That, not the named term, is what this package is for.

Speed is not a reason to prefer one over the other — see Performance.

Performance

gwespmatch maintains a shared-partner cache, so a toggle costs O(degree) rather than O(degree²). This is the same strategy that makes ergm's own gwesp fast, and it puts the two at parity.

The key observation is that a two-path a–k–b survives a nodematch filter exactly when both legs are within-group — i.e. when attr[a] == attr[k] == attr[b], the triad condition. So the triad shared-partner cache is just an ordinary shared-partner cache built on the within-group subgraph, which is precisely what F(~gwesp, ~nodematch) does internally. gwespmatch builds that same cache directly (see src/gwespmatch_spcache.c), with a flag that turns the filter off to get "dyad".

MCMC, undirected, n = 250, 4000 draws

mean degree	gwespmatch	F(~gwesp, ~nodematch)	ratio
15	0.72 s	0.65 s	1.11×
33	0.73 s	0.66 s	1.11×
56	0.90 s	0.84 s	1.07×
89	1.17 s	1.14 s	1.03×

The ratio is flat in density and converges to parity — the two now share a complexity class, and what remains is constant-factor overhead. (Before the cache, the same comparison degraded from 1.12× to 1.56× as degree rose.)

MCMC, directed, n = 200, all five types

type	gwespmatch	F(~dgwesp, ~nodematch)	ratio
OTP	0.75 s	0.60 s	1.25×
ITP	0.55 s	0.48 s	1.15×
RTP	0.39 s	0.37 s	1.05×
OSP	0.38 s	0.39 s	0.97×
ISP	0.39 s	0.38 s	1.03×

homophily = "dyad" runs on the same machinery (the cache is simply built unfiltered), so it gets the same O(degree) behaviour — and there is no F() construction to compare it against.

The summary statistic (summary(nw ~ gwespmatch(...))) is still computed by direct enumeration rather than from the cache. It is called once per model, so this has not been worth optimising; it also serves as an independent check that the cached change statistics are right — the test suite verifies the two agree toggle-by-toggle.

How the cross-group types are cached

The two original caches cannot answer the cross-group types. Both store a single number per dyad — all partners, or partners matching both endpoints — and the cross-group types need to know which endpoint a partner matches, a distinction that single number has already thrown away.

The fix is to bin partners by their own category: cache_c(a, b) = shared partners k of (a, b) with x_k = c (src/gwespmatch_catcache.c). Writing L = cache[x_tail] and R = cache[x_head], every type is then O(1):

type	derivation	caches used
dyad	T	unfiltered
in.in	B	within-group
in.out	T − B	both originals — no new machinery
out.tail / out.head	L / R	per-category
out.one	L + R	per-category
out.none	T − L − R	unfiltered + per-category

Two properties make this the cheap way to do it. First, "k is a shared partner of {a,b} and x_k = c" is symmetric in a and b, so each per-category map keeps exactly the key structure of the corresponding original — the undirected types stay undirected and nothing needs re-keying. (The obvious alternative, filtering by which leg is within-group, is not symmetric and would force UTP/RTP/OSP/ISP onto directed maps.) Second, every two-path lands in exactly one map, so total maintenance per toggle is unchanged from a single unfiltered cache.

dyad and in.in were left untouched and still use the original caches alone.

Measured, on undirected random networks with n = 300 and three attribute categories, timing 5×10⁴ MCMC steps — cached lookup against recounting the same statistic:

mean degree	in.out	out.one	out.none
6	0.60 s vs 2.47 s (4.1×)	0.60 s vs 6.39 s (10.6×)	1.32 s vs 4.17 s (3.2×)
18	1.28 s vs 6.47 s (5.1×)	1.26 s vs 9.42 s (7.5×)	1.42 s vs 4.37 s (3.1×)
53	1.92 s vs 7.84 s (4.1×)	2.36 s vs 14.72 s (6.2×)	2.59 s vs 8.02 s (3.1×)
107	3.00 s vs 17.66 s (5.9×)	3.30 s vs 60.76 s (18.4×)	3.49 s vs 15.72 s (4.5×)

Both builds maintain the caches, so this isolates the lookup rather than flattering it. The widening gap is the point: dropping from O(degree²) to O(degree) is worth ~3× on sparse data and grows from there.

The complementary check is that the cross-group types now track the plain cached path instead of drifting away from it — out.one relative to dyad runs 1.50 → 2.23 → 1.34 → 1.29 across those same four densities. Flat-to-falling is what "same complexity class" looks like.

Note that ergm.godfather is useless for measuring any of this: at ~57 µs per toggle it is dominated by R overhead, which swamps the change statistic entirely. Time simulate() with a large MCMC.burnin instead.

A note on RTP

gwespmatch implements ergm's documented RTP rule (k is an RTP partner of (i,j) iff i ↔ k ↔ j), and computes it directly — it does not use ergm's shared-partner cache.

While building this package we found that ergm's shared-partner cache miscomputed edgewise RTP (wrong on ~18% of random directed networks), which we reported and fixed upstream in statnet/ergm#656.

• ergm versions without that fix: dgwesp(type = "RTP") and desp(type = "RTP") are unreliable by default. They agree with the documentation — and with this package — only when called with term.options = list(cache.sp = FALSE). Use that if you compare the two.
• ergm versions with the fix: no special handling needed; the default path agrees with gwespmatch.

Either way gwespmatch is unaffected: its RTP numbers were always correct.

References

• Agneessens, F., Lomi, A. & Pistocchi, V. (2026). The Role of Third Parties in Homophilous Tendencies: Three Competing Theories. 9th European Conference on Social Networks (EUSN 2026), Norrköping.
• Grund, T. U. & Densley, J. A. (2015). Ethnic Homophily and Triad Closure: Mapping Internal Gang Structure Using Exponential Random Graph Models. Journal of Contemporary Criminal Justice, 31(3), 354–370. doi:10.1177/1043986214553377
• Hong, C.-S., Paik, A., Ballakrishnen, S., Silver, C. & Boutcher, S. (2024). Categorical closure: Transitivity and identities in longitudinal networks. Social Networks, 79, 76–92. doi:10.1016/j.socnet.2024.06.004
• Hunter, D. R. (2007). Curved exponential family models for social networks. Social Networks, 29(2), 216–230. doi:10.1016/j.socnet.2006.08.005

Author

Benjamin Rosche (benrosche@nyu.edu) · GPL-3