Fitting models#

ergmx.ergm() fits a model to a network:

fit = ergmx.ergm(network, formula, seed=1)

How it fits depends on whether the model’s terms are dyad-independent: whether the statistics’ change from adding a tie depends only on that tie (edges, nodematch, nodecov…) or also on the rest of the network (mutual, triangle, gwesp…).

Dyad-independent models: the exact MLE#

Without dependence, an ERGM is a logistic regression on the dyads, and ergmx computes its maximum likelihood estimate exactly, with the log-likelihood, AIC and BIC. Do wealthier Florentine families marry more?

import ergmx
from ergmx import datasets

flomarriage = datasets.load("flomarriage")
fit = ergmx.ergm(flomarriage, "edges + nodecov('wealth')")
fit.summary()

Maximum Likelihood Results:

                 Estimate  Std. Error  MCMC %  z value  Pr(>|z|)
edges             -2.5949      0.5361       0   -4.841    <1e-04 ***
nodecov.wealth     0.0105      0.0047       0    2.256   0.02406 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log-likelihood: -51.5543   AIC: 107.1086   BIC: 112.6836

Dyad-dependent models: the Monte Carlo MLE#

With dependence, the likelihood has a normalizing constant that sums over every possible network, so ergmx estimates it by MCMC, as ergm does:

  1. It starts from the maximum pseudo-likelihood estimate (MPLE).

  2. It simulates networks at the current coefficients, in parallel chains, and moves the coefficients towards values whose simulated networks have, on average, the observed statistics (Hummel et al. 2012 step lengths with a log-normal approximation).

  3. It stops when two consecutive steps are full steps, then draws a larger final sample for the standard errors.

mesa = datasets.load("faux.mesa.high")
fit = ergmx.ergm(
    mesa, "edges + nodematch('Grade') + nodematch('Race') + gwesp(0.5, fixed=TRUE)", seed=1
)
fit.summary()

Monte Carlo Maximum Likelihood Results:

                  Estimate  Std. Error  MCMC %  z value  Pr(>|z|)
edges              -6.3235      0.1541       0  -41.045    <1e-04 ***
nodematch.Grade     1.9798      0.1748       0   11.327    <1e-04 ***
nodematch.Race      0.2658      0.1184       0    2.244   0.02482 *
gwesp.fixed.0.5     1.2354      0.0831       0   14.864    <1e-04 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log-likelihood: -868.4783 (MC SE 0.222)   AIC: 1744.9566   BIC: 1776.7485
Converged after 6 iterations (4 chains, 1024 samples).

The standard errors include the Monte Carlo error of the estimate; the MCMC % column is its share of the variance, as in ergm. A last line says whether the estimation converged. The log-likelihood is estimated too: see Model comparison.

The estimates are also available directly:

fit.coef

{'edges': -6.323542379477252,
 'nodematch.Grade': 1.9797518041863202,
 'nodematch.Race': 0.2657610639159213,
 'gwesp.fixed.0.5': 1.2354403176948652}

fit.stderr

{'edges': 0.1540651294125193,
 'nodematch.Grade': 0.17478930077864277,
 'nodematch.Race': 0.11842283149765728,
 'gwesp.fixed.0.5': 0.08311865912584615}

fit.params and fit.cov give the same as an array and a matrix, in the order of fit.names.

Seeds and threads#

The MCMC is random: pass seed= for reproducible results. Chains run in parallel threads, four by default (or one per CPU, if fewer), and a fit is reproducible for a given seed and number of chains. n_chains=1 runs on a single thread.

To watch the estimation, turn on logging:

import logging

logging.basicConfig(level=logging.INFO, format="%(message)s")
fit = ergmx.ergm(mesa, "edges + nodematch('Grade') + gwesp(0.5, fixed=TRUE)", seed=2)
logging.getLogger().setLevel(logging.WARNING)

iteration 1: interval 1024, effective size 69, step length 0.72
iteration 2: interval 1024, effective size 39, step length 1.00
iteration 3: interval 2048, effective size 42, step length 1.00
iteration 4: interval 4096, effective size 102, step length 1.00
iteration 5: interval 4096, effective size 111, step length 1.00
final iteration: interval 4096, step length 1.00, p-value 0.328

Each iteration reports the MCMC interval (proposals between sampled networks), the effective sample size and the step length. When the samples are too autocorrelated, the interval grows.

Settings#

ergmx.Control holds the settings of the MCMC and the estimation, with ergm’s defaults. Pass a Control, or override single settings as keyword arguments:

ergmx.ergm(mesa, formula, samplesize=2048, n_chains=8)
ergmx.ergm(mesa, formula, control=ergmx.Control(interval=4096))

The most useful are samplesize (networks per iteration), interval, burnin, n_chains and effective_size (the effective sample size the interval adapts to).

Other estimates and starting values#

estimate="MPLE" stops at the MPLE, which is quick but has unreliable standard errors for dyad-dependent models. estimate="CD" stops at the contrastive divergence estimate: the coefficients at which networks a few MCMC steps away from the observed one have, on average, its statistics. It has no standard errors.

init sets where the Monte Carlo MLE starts: "MPLE" (the default), "CD", or your own coefficients, such as those of a previous fit:

ergmx.ergm(mesa, formula, init="CD")
ergmx.ergm(mesa, formula, init=previous_fit.params)

When a fit doesn’t converge, or stops with a DegeneracyError, see Convergence and MCMC diagnostics.

Estimating the decay: curved models#

gwesp(0.5, fixed=TRUE) fixes the decay at 0.5. Without fixed=TRUE, as is ergm’s default, the decay is estimated with the other coefficients, starting from 0.5, and the model is a curved exponential family:

curved = ergmx.ergm(
    mesa, "edges + nodematch('Grade') + nodematch('Race') + gwesp(0.5)", seed=1
)
curved.summary()

Monte Carlo Maximum Likelihood Results:

                  Estimate  Std. Error  MCMC %  z value  Pr(>|z|)
edges              -6.3659      0.1623       0  -39.220    <1e-04 ***
nodematch.Grade     1.9763      0.1736       0   11.381    <1e-04 ***
nodematch.Race      0.2715      0.1164       0    2.332   0.01969 *
gwesp               1.3441      0.1373       0    9.788    <1e-04 ***
gwesp.decay         0.3815      0.1053       0    3.624   0.00029 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log-likelihood: -868.2745 (MC SE 0.192)   AIC: 1746.5491   BIC: 1786.2890
Converged after 6 iterations (4 chains, 1024 samples).

gwesp.decay is the estimated decay, with its standard error. The model’s statistics are then the counts of ties with each number of shared partners (curved.observed), which the two parameters weight: see the term reference. R’s ergm fits this model only from a contrastive divergence start (init.method = "CD"): from its default start, a simulated network exceeds its cutoff of 30 shared partners and the fit stops with an error. ergmx counts those ties in an overflow statistic, and its estimates match R’s.

Fixed coefficients#

offset(term) fixes a term’s coefficients instead of estimating them, at the values given as offset_coef, in formula order:

fixed = ergmx.ergm(flomarriage, "offset(edges) + nodecov('wealth')", offset_coef=[-2.6])
fixed.summary()

Maximum Likelihood Results:

                 Estimate  Std. Error  MCMC %  z value  Pr(>|z|)
offset(edges)     -2.6000                                         (offset)
nodecov.wealth     0.0106      0.0022       0    4.838    <1e-04 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log-likelihood: -51.5543   AIC: 105.1087   BIC: 107.8962

Offsets don’t count as parameters in AIC and BIC. A coefficient of -inf forbids the ties the term counts: see Fixed coefficients.

Constraints and missing ties#

constraints= restricts the networks the model puts probability on, for example to bounded degrees in fixed-choice designs: see Constraints. Dyads whose value is unknown are marked as edges with na=True; the fit is then conditional on the observed ones: see Missing ties.

Saving fits#

A fit takes seconds to minutes; fit.save(path) keeps it, with its network, estimates, MCMC sample and settings, and ergmx.load_fit() brings it back, to summarize, simulate, check or compare as before (R’s saveRDS() and readRDS()):

fit.save("mesa_fit.pkl")
fit = ergmx.load_fit("mesa_fit.pkl")

The file is a Python pickle (the fits pickle, too, with pickle.dump()), so load only files you trust; loading a file saved by another version of ergmx warns.