# Validation

`ergmx` is tested against R's ergm 4.12, ergm.multi 0.3.0 and tergm 4.2.2 on
ergm's own networks (flomarriage, samplk1 to samplk3, faux.mesa.high and
faux.dixon.high, the latter two also with missing dyads), multinets'
`linked_sim`, Davis's Southern Women, a simulated bipartite network,
Sampson's monks as a sample and as a series of networks, and ergm.multi's
225 weekday household networks (Goeyvaerts), and against exact results on
networks small enough to enumerate every possible network. The R scripts in
`scripts/` store R's results in `tests/data/`, and the test suite compares.

## Against R's ergm

| Check | Result |
|---|---|
| Statistics of ergm's terms and operators, with their options (`levels=`, `by=`, `homophily=`, `attr=`, `nodes=`...), and interactions, 97 models | identical to R's `summary()` (to 1e-12), names included, except `transitive`, `intransitive`, `dyadcov`'s `utri` and `ltri`, and the edgewise RTP statistics (below) |
| MPLE, 60 models, with offsets, interactions, `F()`, `S()`, `N()` (with `subset`, `offset` and `label`), tergm's operators (also with missing dyads imputed by `NA.impute`), `blocks`, `Dyads`, `fixedas`, `fixallbut`, `blockdiag`, bipartite networks and curved terms | identical to R (to 1e-6; curved models to 1e-3, where R's optimizer stops on a flat optimum: ergmx's pseudo-likelihood is at least R's) |
| Dyad-independent MLE, standard errors, log-likelihood and BIC, 25 models, with offsets, interactions, `N()`'s `subset` and `offset`, `NA.impute`, `blocks`, `Dyads`, `fixedas`, `fixallbut`, `blockdiag`, `S()`, missing dyads, bipartite networks, samples and series of networks | identical to R (to 1e-6; standard errors to 1e-3, the tolerance of R's `glm`) |
| Monte Carlo MLE, 7 models (4 undirected, 3 directed), 3 to 10 seeds each | within 0.12 standard errors of R's estimates; standard errors 0.89 to 1.12 times R's |
| Monte Carlo MLE with constraints (`bd`, `blocks`, `degrees`, `odegrees`), missing dyads, offsets, `F()`, `esp` and multilevel models, 12 models x 5 seeds | all converged, within 0.17 standard errors of R's estimates; standard errors 0.91 to 1.09 times R's |
| Monte Carlo MLE of `concurrent`, `twopath`, OSP shared partners, bipartite models and curved models (gwesp, directed and undirected, and gwb1degree), 10 models x 3 to 5 seeds | all converged, within 0.2 standard errors of R's estimates; standard errors 0.90 to 1.12 times R's, except the curved bipartite model, whose likelihood is nearly flat in the decay (R's standard error of the decay is twice its estimate) |
| Monte Carlo MLE of samples of networks (Sampson's monks; 225 households with `N()` linear models) and of series (tergm's CMLE, one and two transitions), 4 models x 5 seeds | all converged, within 0.10 standard errors of R's estimates; standard errors 0.94 to 1.05 times R's |
| Monte Carlo MLE with `N()`'s `subset` and `label`, and of a series with missing dyads imputed (`NA.impute="next"`) and missing in the networks transitioned to, 2 models x 5 seeds | all converged, within 0.09 standard errors of R's estimates; standard errors 0.93 to 1.03 times R's |
| Monte Carlo MLE of a multilevel model with `S()` (each level and the ties between them, with gwesp and gwb1dsp), 5 seeds | all converged, within 0.08 standard errors of R's estimates; standard errors 0.96 to 1.04 times R's |
| Monte Carlo MLE with `degree(by=)`, and the constraints `edges`, `b1degrees` and `bd(attribs=)`, 4 models x 5 seeds | all converged, within 0.16 standard errors of R's estimates; standard errors 0.93 to 1.11 times R's |
| Goodness of fit, directed and undirected | observed distributions and p-values identical to R's; simulated distributions agree within Monte Carlo error |
| ergm.multi's `gofN()`, 225 households, the model's statistics and five others, 2,000 simulations | observed statistics identical to R's (but `degree0` and `isolates`, below); fitted values, variances and Pearson residuals agree within R's Monte Carlo error |
| tergm's EGMME, formation and persistence of edges with a mean duration, and with `degree(1)` too | within 0.4 standard errors of the mean of R's estimates over 3 seeds; standard errors within R's range, which spans a factor of 1.7 between its seeds; the edges model's estimate is also the exact one (below) |
| Conditional tie probabilities (`predict()`), 5 models: undirected, directed, missing dyads, curved, 20 household networks | identical to R's `predict()` (to 1e-12), on the same dyads (R's formula method leaves out missing dyads) |
| Average marginal effects, 4 models | identical to ergMargins' (to its 5 significant digits); the standard errors match a numerical delta method (to 1e-5), and ergMargins' when the probabilities are held fixed, as it holds them |
| Confidence intervals, 4 models | identical to R's `confint()` (to 1e-3, the accuracy of R's `glm` standard errors) |
| Tables of results | identical to texreg's `screenreg()`, `texreg()` and `htmlreg()`, character for character, on exactly fitted models |
| Contrastive divergence | a fixed point of its defining equation; the Monte Carlo MLE from a CD start matches R's |

R's own standard errors vary by about 15% between seeds on the small
networks, so the comparison of standard errors is only as tight as R allows.
ergm imputes missing dyads at random before its MPLE, so MPLEs of
dyad-dependent models with missing dyads are not compared. R fits two of the
curved models only from a contrastive divergence start; ergmx fits them from
its default start. ergm.multi and tergm can't fit curved terms inside `N()`
or `Form()` (their MPLE's gradient is not finite, and contrastive divergence
stops with an error), so only their statistics are compared; ergmx's fits of
them are checked below.

### Discrepancies in ergm 4.12.0

- **Edgewise RTP statistics.** ergm 4.12.0's `esp`, `gwesp` and `nsp` with
  `type = "RTP"` depend on the order of the vertices: relabeling
  faux.dixon.high changes `summary(net ~ esp(0:3, type = "RTP"))` from 889,
  269, 37, 2 to 886, 273, 36, 2. Its shared-partner cache, on by default,
  reads them with the wrong key; ergm's development version fixes it
  ([statnet/ergm#656](https://github.com/statnet/ergm/pull/656)). With the
  cache off (`term.options = list(cache.sp = FALSE)`), ergm 4.12.0 gives 847,
  309, 39, 2, as the definition computed in R does, and ergmx's `esp`,
  `gwesp`, `nsp` and `gwnsp` of type RTP are identical to ergm's.
- **`transitive`.** ergm documents it as the number of transitive triads
  (types 030T, 120D, 120U and 300: 371 on faux.dixon.high, 13 on samplk3,
  counted by R's igraph) but computes the number of transitive triples, the
  same as `ttriple` (1,254 and 49). ergmx counts the triads, which match the
  triad census on every directed network of 4 vertices, and warns that ergm
  differs; `ttriple` reproduces ergm's statistic.
- **`intransitive`.** The same: documented as the number of intransitive
  triads (types 111D, 201, 111U, 021C and 030C: 5,991 on faux.dixon.high,
  70 on samplk3, by igraph's triad census) but computed as the number of
  intransitive triples, twopath minus ttriple (6,752 and 92). ergmx counts
  the triads, which match R's `triadcensus` and igraph's, and warns.
- **`dyadcov`.** Documented as the covariate summed over mutual,
  upper-triangular asymmetric and lower-triangular asymmetric dyads; on a
  network whose only tie goes from vertex 1 to vertex 2, in the upper
  triangle, ergm's `utri` is 0 and its `ltri` 1. ergmx follows the
  documentation, so its `utri` and `ltri` are ergm's `ltri` and `utri`, and
  warns.
- **`gofN()`'s `degree0` and `isolates`.** ergm.multi's `gofN()` leaves out
  the statistics of the empty network, which for `degree0` and `isolates` are
  the network size: in the households, it reports the observed and fitted
  isolates as their counts minus the size (-4 for a household of four
  without isolates). The Pearson residuals are unaffected. ergmx reports the
  counts, and warns.
- **`smalldiff`** is documented as counting differences less than the
  cutoff, but its code counts those at most the cutoff (193 ties of
  faux.mesa.high with grades at most 2 apart, against 178 strictly less);
  ergmx follows the code, as the argument's description, "maximum", does.

## Against exact results

| Check | Result |
|---|---|
| MCMC stationary distribution, with 3 mixtures of TNT and triadic proposals | expected statistics match exact enumeration of every network on 3 and 6 vertices (undirected) and 4 vertices (directed) |
| Constrained MCMC: `bd` (directed and undirected, also by alter class), `degrees` (both), `odegrees`, `idegrees`, `b1degrees`, `b2degrees`, `edges`, `blocks`, and sampling conditional on observed dyads | expected statistics match exact enumeration of every network each constraint allows; a test that only reversing cyclic triples can pass checks that directed `degrees` reaches every network |
| `-inf` offsets | the MLE equals the logistic regression without the forbidden dyads |
| Bipartite MCMC | expected statistics match exact enumeration of the 512 networks of a 3 x 3 bipartite network |
| Curved terms | eta . counts equals theta times the fixed-decay statistic exactly, and the Jacobian matches finite differences; inside `N()`, with a linear model, too |
| Several networks: `Networks()` and `NetSeries()` | the MCMC matches exact enumeration of every network with ties within the networks (two undirected networks; two transitions of directed networks, with discordant-dyad and triadic proposals), and never ties networks together |
| Curved terms inside `N()` and `Form()` | at the curved MPLE's decay, the fixed-decay MPLE has the same coefficients (to 1e-8); a series simulated with decay 0.7 gives back 0.64 |
| Dynamic simulation | one time step is a draw from the transition's model (exact enumeration); with tergm's stopping rule, ties form and dissolve at the dyad-independent model's exact rates (within 5% and 8%) |
| Dyad-independent CMLE | the log-odds that non-ties became ties and that ties persisted, exactly |
| Forward simulation with time trends (`lm=~.Time`) | each step's formation and persistence rates are those predicted for its time, within Monte Carlo error |
| `NA.impute` | each option fills the dyads it can, in turn, as tergm's (`next` from the next wave, also through a wave missing them too) |
| `N()`'s `offset` | the offset statistics' coefficients are 1 and the Jacobian matches finite differences, for curved terms too; the MCMC with `subset` and `offset` matches exact enumeration |
| `gofN()` | each network's fitted value and variance of its edges are the dyad-independent model's exact ones (within Monte Carlo error); with missing dyads, the observed value and its variance are the exact conditional ones |
| EGMME | the edges model, whose equilibrium is known (formation and persistence probabilities that give the density and mean duration), within 0.5 standard errors; its standard errors are the delta method's from the exact stationary covariances; tie ages follow every tie of a simulation |
| The sampler's tracked statistics for every new term (degree ranges, by attribute and with homophily, triad census, trails, Simmelian and transitive ties, covariate ranges, distinct neighbour types, the bipartite terms...) | equal the statistics recomputed from scratch on the sampled networks (to 1e-9), so the change statistics of removing ties are right too |
| Constraints `edges`, `b1degrees`, `b2degrees`, `bd(attribs=)`, `fixedas`, `fixallbut`, `Dyads`, `blockdiag`, `observed` | every simulated network keeps what the constraint fixes, and moves elsewhere |
| MPNet's directed configurations (35, from its manual) and EXTA, EXTB and ASAXASB | equal to their definitions computed from adjacency matrices (to 1e-12) on random directed two-level networks; the MCMC matches exact enumeration of every directed network on 2 + 2 vertices with affiliations from A to B; `S()` between two sets of a directed network matches R's |
| Estimated decays of the multilevel configurations (18 terms) | eta . counts equals theta times the fixed-decay statistic exactly; at the curved MPLE's decay, the fixed-decay MPLE has the same coefficients; networks simulated with a decay of 0.7 give it back within 1.1 standard errors |
| Goodness of fit by level | the observed distributions within each level and of the affiliations equal igraph's |
| MPNet's 16 multilevel configurations, which no R package has | equal to their definitions computed from the levels' adjacency matrices (to 1e-12), on random two-level networks and `linked_sim`, for three decays; the MCMC matches exact enumeration of every network on 6 vertices of two levels (and one of neither), with S() too |
| Log-likelihood, directed and undirected | unbiased against exact enumeration (6 and 4 vertices), with standard errors that match the spread across seeds |

## The log-likelihood

For 5 dyad-dependent models, `ergmx`'s log-likelihood is within one standard
error of a high-precision estimate (128 points along the path, 4 times the
samples). R's ergm, which integrates with a 16-point midpoint rule, is off by
0.1 to 2.1 units:

| Model | R ergm | High-precision | ergmx (3 seeds) |
|---|---|---|---|
| samplk3, edges + mutual | -133.99 | -133.88 | -133.86, -133.90, -133.87 |
| samplk3, + gwesp(0.5) | -131.36 | -131.04 | -130.99, -131.04, -131.15 |
| faux.mesa.high, gwesp(0.5) | -867.56 | -867.09 | -867.13, -867.15, -867.09 |
| faux.mesa.high, gwdegree + gwesp | -866.01 | -865.27 | -865.36, -865.20, -865.74 |
| faux.dixon.high (directed), gwesp(0.1) | -4206.10 | -4208.20 | -4208.56, -4208.57, -4207.76 |

## Performance

`benchmarks/benchmark.R` and `benchmarks/benchmark.py` fit three models with
both packages' defaults, which include the log-likelihood; median of 3 seeds
on an Apple M4 Pro. R's ergm runs on one thread, and the single-threaded
ergmx run is limited to one thread too.

| Model | Vertices | R ergm 4.12 | ergmx, 1 thread | ergmx, all threads | Largest difference |
|---|---|---|---|---|---|
| faux.mesa.high, gwesp(0.5) | 205 | 16.8 s | 10.1 s (1.7x) | 2.4 s (7.1x) | 0.04 SE |
| faux.magnolia.high, gwesp(0.25) | 1,461 | 17.3 s | 14.0 s (1.2x) | 3.3 s (5.2x) | 0.07 SE |
| faux.dixon.high (directed), gwesp(0.1) | 248 | 142 s | 80 s (1.8x) | 17 s (8.6x) | 0.07 SE |
| 225 household networks (ergm.multi), with `N()` linear models | 2 to 7 each | 23.3 s | 13.2 s (1.8x) | 2.9 s (8.0x) | 0.06 SE |

The last column is the largest difference between the two packages'
estimates, in R's standard errors (the households: one R run, from
`scripts/r_reference.R`). On small models, such as the series of Sampson's
monks, both take about a second. `ergmx`'s log-likelihood samples about 17
times more than ergm's, which is what makes it accurate; without the
log-likelihood on either side, a single thread was 1.9 to 3.6 times faster
than R on these models.

### Larger networks

`benchmarks/make_scale.py`, `scale.R` and `scale.py` fit the magnolia model
to faux.magnolia.high and to networks like it of 2,500 to 10,000 vertices
(simulated from R's fit, with the same mean degree), and a model with
homophily and gwesp to 500 classrooms of 20 students, combined with
`Networks()` (`N()` in ergm.multi); one seed each, the same machine:

| Network | Vertices | MPLE: R | ergmx | MLE: R | ergmx, 1 thread | ergmx, all threads | Largest difference |
|---|---|---|---|---|---|---|---|
| faux.magnolia.high | 1,461 | 0.2 s | 0.01 s | 17.2 s | 13.9 s (1.2x) | 3.2 s (5.5x) | 0.10 SE |
| magnolia-like | 2,500 | 0.4 s | 0.02 s | 11.4 s | 6.6 s (1.7x) | 1.4 s (8.4x) | 0.06 SE |
| magnolia-like | 5,000 | 1.4 s | 0.08 s | 12.0 s | 26.8 s (0.4x) | 5.8 s (2.1x) | 0.11 SE |
| magnolia-like | 10,000 | 5.6 s | 0.3 s | 36.9 s | 27.5 s (1.3x) | 5.6 s (6.6x) | 0.07 SE |
| 500 classrooms | 10,000 | 2.6 s | 0.02 s | 91.3 s | 59.7 s (1.5x) | 15.3 s (6.0x) | 0.02 SE |

The MPLE is 17 to 110 times faster: its data are the distinct rows of change
statistics, built in parallel threads, and the logistic regression runs on
those. The Monte Carlo MLEs take a varying number of iterations: on 5,000
vertices, ergmx's took 8 where the others took 4 or 5, and half of each
fit's time is the log-likelihood. Memory stays below 0.6 GB for all of them.