Exact Likelihood Inference for Snowball-Sampled Erdős-Rényi Networks

Abstract: Network data obtained through link-tracing designs, such as snowball sampling, are collected through a mechanism that depends on the very structure the analysis seeks to estimate. Ignoring this dependence and treating the observed sample as though it were itself a complete network can lead to substantially biased inference. While the resulting selection problem is intractable in general, we show that it admits an exact solution for $r$-wave snowball samples, with full-neighbourhood recruitment, drawn from an Erdős--Rényi population. We derive the exact likelihood of such a sample and show that it defines a curved exponential family in the edge probability $π$, with a low-dimensional sufficient statistic. Building on this result, we obtain the maximum likelihood estimator of $π$ that correctly accounts for the sampling design and, as a function of the minimal sufficient statistic, makes full use of the information in the sample. Simulation studies show that this correction substantially reduces bias relative to the naive estimator, remaining effectively unbiased even when the sample covers as little as 0.1\% of the network. We further construct valid confidence intervals for $π$ by inverting a test built on the exact sampling distribution, approximated via Monte Carlo simulation. Simulation studies confirm that these confidence intervals attain the nominal coverage level within Monte Carlo error across a range of edge probabilities and numbers of waves.

• The paper derives a closed-form likelihood for r-wave snowball samples from Erdős–Rényi networks, showing that the edge count and a correction term for unsampled vertices form a minimal sufficient statistic.
• The snowball-corrected maximum likelihood estimator, |E^(r)|/T(n), is approximately unbiased even when samples cover fewer than 1% of the population, while naive estimators can overstate sparse network density by one to two orders of magnitude.
• The paper constructs conservative confidence sets by inverting exact-distribution Monte Carlo tests, achieving empirical coverage of 0.943–0.962 and intervals about 44.6% narrower than initial likelihood brackets.

Motivation and problem

Link-tracing designs such as snowball sampling collect network data through the very edges the analysis seeks to estimate, so the sampling mechanism is informative about the population model. Treating a snowball sample as if it were a complete graph induces systematic upward bias in density estimates, since high-degree vertices are over-represented. The general likelihood-based remedy—summing over all completions of the unsampled portion of the network, as formalized by Handcock and Gile's notion of "amenability"—is combinatorially intractable. This paper resolves the base case posed as Research Problem 3.1 by Crane (2018): an exact, closed-form distribution for $r$-wave snowball samples with full-neighbourhood recruitment from an Erdős–Rényi (ER) population.

Exact likelihood and sufficient statistic

The key structural observation is that the joint probability of the sampled adjacency matrix $y^{(r)}$ and the wave sets factorizes into (i) a conditional probability of $y^{(r)}$ given the wave sets, which is constrained by two design rules—no edges between non-adjacent waves, and at least one edge from each wave-$k$ vertex to wave $k-1$—and (ii) a marginal probability of the wave sets, where each wave constitutes a Bernoulli sampling design with success probability $1-(1-\pi)^{n_{k-1}}$ over the remaining vertices. The design-induced factors cancel exactly, yielding

$P(y^{(r)}, V^{(1)},\dots,V^{(r)} \mid V^{(0)}) = \pi^{|E^{(r)}|}(1-\pi)^{T(n)-|E^{(r)}|}\,\mathds{1}\{y^{(r)} \in \mathcal{S}(n)\},$

where $T(n) = |W| + |A| + |J| + U\sum_{k=0}^{r-1} n_k$ counts all constrained vertex pairs: within-wave pairs, adjacent-wave pairs, non-adjacent-wave pairs, and pairs between unsampled vertices ($U$ of them) and waves $0,\dots,r-1$. The last term encodes the informativeness of exclusion—an unsampled vertex must have no tie to any earlier-wave vertex—and this is precisely what the naive estimator discards.

The distribution is a curved exponential family with natural statistic $y^{(r)}$0, which the authors prove is minimal sufficient for $y^{(r)}$1: individual wave sizes are sufficient but not minimal, since distinct wave-size vectors with equal $y^{(r)}$2 carry identical information. The derivation is verified by showing that per-wave simulation and direct snowball sampling from a fully generated ER graph produce indistinguishable joint distributions of $y^{(r)}$3; a permutation test based on the two-sample energy statistic gives a Monte Carlo $y^{(r)}$4-value of 0.678 over 10,000 replicates at $y^{(r)}$5, $y^{(r)}$6.

Snowball-corrected maximum likelihood estimation

The MLE takes the simple form $y^{(r)}$7, differing from the naive estimator only in the denominator: the corrected version replaces the count of non-adjacent-wave pairs inside the sample with the expected number of unsampled pairs. The bias of the naive estimator vanishes only as $y^{(r)}$8.

Simulations with $y^{(r)}$9 and sparse edge probabilities $y^{(r)}$0 show that for one-wave samples the naive estimator exceeds the truth by one to two orders of magnitude, and remains biased even at three waves unless the sample covers more than roughly two-thirds of the population. The corrected estimator is approximately unbiased across all settings, including samples averaging fewer than 100 vertices—about 0.7% of the network. Because $y^{(r)}$1 is a function of the minimal sufficient statistic, it uses all information in the sample; the practical implication is that valid point estimation of edge probability requires neither dense coverage nor large numbers of waves under this model.

Monte Carlo confidence intervals via test inversion

The paper deliberately avoids asymptotic normality arguments, showing that the limiting behavior of $y^{(r)}$2 depends delicately on the relationship between $y^{(r)}$3 and $y^{(r)}$4: with $y^{(r)}$5 fixed and $y^{(r)}$6, wave 2 absorbs nearly the whole population; under $y^{(r)}$7, the first wave converges to a Poisson rather than concentrating. Instead, confidence sets are obtained by inverting a two-sided test built on the exact sampling distribution of $y^{(r)}$8, using conservative importance-sampling p-values (Harrison's add-one correction) evaluated on batches simulated at a proposal value $y^{(r)}$9 and reweighted across candidate values, with redraws triggered when the effective sample size falls below a floor. The resulting set is guaranteed conservative for any $k$0, $k$1, and $k$2; its convex hull is reported, avoiding any quasiconcavity assumption on the p-value function, and initialized from a relative-likelihood bracket located by bisection.

Across 9,000 replicates spanning nine $k$3 combinations, empirical coverage ranges from 0.943 to 0.962 (mean 0.952), consistent with nominal coverage within Monte Carlo error. Interval width shrinks by roughly an order of magnitude per additional wave, and the refined intervals are on average 44.6% narrower than the starting likelihood bracket while retaining the coverage guarantee.

Limitations and open questions

Three restrictions are acknowledged explicitly. First, the framework applies only to the ER model; extension to dyad-independent models conditional on latent structure (latent space, stochastic block, or graphon models) is conjectured to be straightforward but not carried out. Second, all results assume an exogenous ego: selection proportional to degree would inflate $k$4 and bias $k$5 upward, and a closed-form correction appears infeasible because the required normalization involves total degree in the unsampled portion of the graph. Third, the convex hull of the Monte Carlo confidence set may be conservative if the p-value function is not quasiconcave; whether quasiconcavity holds in this setting is asserted but not proven. A useful robustness result is established: unobserved edges within the last wave are handled by simply dropping those $k$6 pairs from numerator and denominator, yielding the adjusted estimator $k$7.

Conclusion

The paper shows that for dyad-independent population models, the snowball sampling mechanism can be incorporated exactly into the likelihood through a low-dimensional correction term accounting for unsampled vertices, resolving Crane's open problem in the ER case. The resulting estimator removes severe finite-sample bias, and exact-distribution-based Monte Carlo intervals achieve nominal coverage without asymptotic approximations. Whether the same closed-form amenability extends to conditionally independent dyad models with latent heterogeneity remains the central open question.