Copula Approximate Bayesian Computation

Using Distribution Random Forests

George Karabatsos111Departments of Educational Psychology, and Mathematics,
Statistics, and Computer Sciences, 1040 W. Harrison St. (MC 147), Chicago,
Illinois 60607, U.S.A. email<EMAIL_ADDRESS>
February 23, 2024

Abstract. This paper introduces a novel Approximate Bayesian Computation (ABC)
framework for estimating the posterior distribution and the maximum likelihood
estimate (MLE) of the parameters of models defined by intractable likelihood
functions. This framework can describe the possibly skewed and high
dimensional posterior distribution by a novel multivariate copula-based
distribution, based on univariate marginal posterior distributions which can
account for skewness and be accurately estimated by Distribution Random
Forests (DRF) while performing automatic summary statistics (covariates)
selection, and on robustly estimated copula dependence parameters. The
framework employs a novel multivariate mode estimator to perform for MLE and
posterior mode estimation, and provides an optional step to perform model
selection from a given set of models with posterior probabilities estimated by
DRF. The posterior distribution estimation accuracy of the ABC framework is
illustrated through simulation studies involving models with analytically
computable posterior distributions, and involving exponential random graph and
mechanistic network models which are each defined by an intractable likelihood
from which it is costly to simulate large network datasets. Also, the
framework is illustrated through analyses of large real-life networks of sizes
ranging between 28,000 to 65.6 million nodes (between 3 million to 1.8 billion
edges), including a large multilayer network with weighted directed edges.

Keywords: Bayesian analysis $\cdot$ Maximum Likelihood $\cdot$ Intractable
likelihood.

##  1 Introduction

Statistical models defined by intractable likelihood functions are important
for the analysis of complex and large datasets which are becoming increasingly
available in scientific fields. The likelihood for a given model may be
intractable either because it does not have an explicit form, the model is
high dimensional (large number of parameters), or because of the large size of
the dataset being analyzed. Approximate Bayesian Computation (ABC) provides
alternative methods for estimating the approximate posterior distribution or
MLE for any model defined by an intractable-likelihood.

Each ABC method is uniquely defined by a specific surrogate approximate model
likelihood, which measures the closeness of (preferably sufficient) summary
statistics of the observed dataset to the same summaries of a random dataset
simulated from the exact likelihood conditionally on any proposed parameters
of the given model under consideration. Indeed, the selection of summary
statistics is an important issue when they are not sufficient for the given
model. The ABC field has received many reviews due to its wide theoretical
scope and applicability (e.g., Sunnåker et al. 2013; Karabatsos & Leisen 2018;
Sisson et al. 2018; Grazian & Fan, 2020; Cranmer et al. 2020; Craiu & Levi,
2022; Karabatsos 2023; Pesonen et al. 2023; Martin et al. 2023).

This paper introduces and studies a novel ABC framework for estimating the
posterior distribution and the maximum likelihood estimate (MLE) of the
parameters of models defined by intractable likelihood functions, which
unifies and extends copula ABC (Li et al. 2017; Nott et al. 2018; Chen &
Gutmann, 2019; Klein et al. 2024), ABC-Random Forests (abcrf) (Raynal et al.
2018), and ABC and other simulation-based (AMLE) methods for calculating the
approximate MLE for intractable likelihood models (Rubio & Johansen, 2013;
Kajihara et al. 2018; Picchini & Anderson, 2017; Gutmann & Corander 2016;
Yildirim et al. 2015; Dean et al. 2014; Gourieroux et al. 1993; McFadden,
1989).

Specifically, the new ABC framework combines: copula ABC’s ability to
approximate high-dimensional posterior distributions based on estimating one-
dimensional marginal posterior distributions (unlike abcrf and AMLE) while
separately estimating the dependence structure of the full multivariate
posterior; with abcrf’s ability to estimate these one-dimensional marginals
(dependent variables) using Distribution Random Forests (DRF) (Ćevid et al.
2022) conditionally on observed data summaries while being able to
automatically select summary statistics (covariates) that are relevant for the
given model under consideration (unlike nearly all copula ABC methods, and
AMLE); and with the ability to calculate the approximate MLE (unlike abcrf and
copula ABC).

Using the ABC framework, the aim is to estimate the posterior distribution of
possibly high-dimensional model parameters at least approximately-well, while
accounting for any skewness in the posterior, and for the fact that the
posterior converges to a multivariate normal distribution as the data sample
size increases, even for misspecified models, under regularity conditions
(Kleijn & van der Vaart, 2012). Also, for the purposes of estimating the
posterior mode and MLE, the (approximate) posterior distribution estimate
should be easy to sample from and compute its probability density function
(pdf). With these requirements, for any given model under consideration, we
approximate the posterior distribution of $d$-dimensional parameters $\theta$
using a tractable multivariate copula-based distribution, estimated from a
Reference Table of samples from the prior distribution of $\theta$ and
corresponding summary statistics of datasets simulated from the exact
likelihood of the given model under consideration.

According to Sklar’s (1959) theorem, a continuous $d$-variate random parameter
vector $\theta$ can be described by a joint (posterior) cumulative
distribution function (cdf), called a copula, which can be uniquely
represented by its univariate marginal densities (distributions) and a copula
density function describing its dependence structure. If the random vector
$\theta$ has any discrete-valued parameters (integer-valued without loss of
generality), then these discrete parameters can be continued by jittering
(Denuit et al. 2005; Madsen & Fang, 2011) in order to maintain a unique
representation of the joint copula-based posterior cdf of $\theta$, while
ensuring that no information is lost and preserving the same dependence
relationship among the model parameters.

Approximating a possibly skewed high-dimensional posterior distribution
requires a suitable choice of multivariate copula-based distribution family.
Choices are limited because it is difficult to construct high-dimensional
copulas (Nelsen et al. 2006, p.105).

Meanwhile, skewness can be introduced into a multivariate copula-based
distribution in either of the following three ways: (1) using skewed
univariate margins (e.g., Hutson et al. 2015; Smith & Vahey, 2016; Ballien et
al. 2022); (2) using a skewed copula density function to account for
asymmetric dependencies in the multivariate distribution, such as a skewed,
grouped Student-$t$, or mixture of $K$ Gaussian or $t$ copula density
functions (e.g.,Wei et al. 2016; Yoshiba 2018; Demarta & McNeil, 2005; Daul et
al. 2003; Kosmidis & Karlis, 2016), or using a covariate-dependent copula
density function (Smith & Klein, 2021; Klein et al. 2024; Acar et al. 2013);
or (3) by combining both skewed marginals and skewed densities, but this can
unnecessarily increase the complexity of the copula model (Ballien et al.
2022). In high dimensions, estimating the parameters of skewed copulas is
challenging and cumbersome (Hintz et al. 2022b), and arguably the mixture
copulas are over-parameterized, with $Kd(d-1)/2$ correlation matrix
parameters, and the number of correlation parameters grows quadratically with
$d$ while it is not easy to reliably estimate $K$ correlation matrices
(Dellaportas & Tsionas, 2019; Qu & Lu, 2021). Besides, the ABC framework aims
to estimate a multivariate posterior mode, and the MLE as a mode of a
posterior-to-prior pdf ratio, while the mode is robust to the shapes of the
tails of the distribution.

Therefore, all things considered, the proposed ABC framework employs a novel
multivariate meta-$t$ (Fang et al. 2002) copula-based distribution
approximation of the posterior, which allows for multivariate skewness by
using skewed marginals which are modeled flexibly and accurately estimated by
DRFs, respectively. This multivariate copula-based distribution is defined by
marginals which are (covariate-dependent) regressions and a multivariate $t$
parametric copula that is independent of covariates, extending previous
related copula models (Pitt et al. 2006; Song, 2000) through the use of DRF.
The dependence structure of the meta-$t$ distribution is defined by the
density of a multivariate $t$ distribution with (posterior) degrees of freedom
and correlation matrix parameters, which allows for symmetric tail dependence
and has robust properties of the $t$ distribution, which inversely-weights
each observation according to its Mahalanobis distance from the center of the
data distribution (Lange et al. 1989).

DRF more appropriately targets the whole distribution of the dependent
variable, instead of mean functions of the dependent variable as done by
previous most earlier RF algorithms and abcrf. The capability of DRF to
automatically select summary statistics enables the new ABC framework to avoid
the potential issues involved with pre-selecting summary statistics (Drovandi
et al. 2023), and to avoid having to introduce an extra prior distribution to
do this selection as done in other ABC methods. Further, for scenarios where
multiple models are considered to analyze the same dataset, the new ABC
framework employs an optional step which can be used to perform model
selection based on posterior model probabilities that are estimated using DRF
(conditionally on the observed data summaries), which provides a simple
extension of previous RF approach to model selection in ABC (Pudlo et al.
2016; Marin et al. 2018) by directly targeting the distribution of the model
index dependent variable, instead of running multiple RFs for binary indices
of models.

Finally, the ABC framework calculates an estimate of the multivariate
posterior mode and MLE of the parameters $\theta$ of the (selected) model. The
MLE is the multivariate mode of the ratio of the meta-$t$ posterior pdf
estimate to the prior pdf, or the posterior mode under a proper uniform prior
supporting the MLE. The new multivariate mode estimator maximizes a meta-$t$
posterior density estimate (or finds the MLE by maximizing the ratio of the
posterior density estimate to the prior density) over parameter samples of
points from the Reference Table (or over an extra set of parameter samples of
points drawn from the meta-$t$ posterior density estimate if necessary), which
avoids a possibly costly grid search involved in multivariate mode estimation
by global maximization. Similarly, other ‘indirect’ multivariate mode
estimation methods (Devroye 1979; Abraham et al. 2003 Abraham et al. 2004; Hsu
& Wu, 2013; Dasgupta & Kpotufe, 2014) use the data sample points to estimate
the mode from a nonparametric kernel, $k$-nearest-neighbor, or parametric
density estimate, while respectively, relying on choice of bandwidth, or $k$,
or on a successful normal transformation of the data. For posterior mode or
MLE estimation, a novel semiparametric mode estimator is being used here,
novel in that it is copula-based and thus applicable to high dimensions as
shown by copulaABC research (Nott et al. 2018), while other multivariate mode
estimators seem only applicable to lower dimensions (see Chacón, 2020, for a
review). The new ABC framework estimates the MLE (mode) while performing
automatic summary statistics (variable) selection via DRFs for the meta-$t$
marginals.

Next, Section §2 describes the proposed ABC framework, summarized in Algorithm
1 below. §3 illustrates the viability of this framework through simulation
studies. One simulation study considers models with analytically computable
posterior distributions, to allow their direct comparisons with posterior
distributions estimated in the ABC framework. The subsequent simulation
studies and real data analyses are based on models defined by intractable
likelihoods. Most models defined by intractable likelihoods have the two
properties that (1) it is possible to rapidly simulate datasets from the
model’s likelihood, and (2) the MLE or another point estimator and its
sampling distribution can be computed for the model. In case (1) ABC methods
and synthetic likelihood methods (Picchini et al. 2023, and references
therein) can be used to estimate the (approximate) posterior distribution and
MLE of the model parameters. In case (2), it is possible to estimate its
approximate posterior distribution of the model parameters using certain ABC
or bootstrap methods which rely on the availability of point estimators (e.g.,
Karabatsos 2023; Newton et al. 2021; Barrientos & Peña 2020; Lyddon et al.
2019; Zhu et al. 2016).

Therefore, this paper focuses on models defined by intractable likelihoods
which do not possess these two properties, at least in large data scenarios.
This is the case for models for network datasets which are ubiquitous in
scientific fields, namely, Exponential Random Graph Models (ERGMs) and
mechanistic network models for large network data settings, where each model
is defined by an intractable likelihood, and it is costly or prohibitive to
simulate large network datasets from the likelihood and to compute the MLE or
other point estimator of the model parameters from large networks. Using these
models, the ABC framework is illustrated through analyses of three large real-
life networks of sizes ranging between twenty thousand to 6.6 million nodes,
including a large multilayer network with weighted directed edges. It can be
concluded that the proposed ABC framework provides practical statistical
methods for analyzing complex and large datasets which are becoming more
widely available.

##  2 ABC Framework

The proposed ABC framework can be summarized in Algorithm 1, shown below. The
algorithm can be applied to any set of $M\geq 1$ models. It outputs an
estimate the ABC posterior distribution
$\pi_{\widehat{m}}^{\text{ABC}}(\theta\mid s(x))$, a probability density
function (pdf) based on corresponding cumulative distribution function (cdf)
$\Pi_{\widehat{m}}^{\text{ABC}}$, the selected best predicted model,
$\widehat{m}\in\\{m\\}_{m=1}^{M}$ among $M\geq 1$ models compared, and the
model’s prior distribution (pdf)$\ \pi_{\widehat{m}}(\theta)$ and its
approximate likelihood of the generic form
${\textstyle\int}f_{\widehat{m},s}(s(x)+\varepsilon\nu\mid\theta)K(\nu)\mathrm{d}\nu$
based on summary statistics vector $s$ and some kernel density function
$\mathcal{K}$.

Step 2, used if multiple models ($M>1$) are considered by the data analyst,
performs model selection using DRF (Ćevid et al. 2022), which extends previous
a Random Forest approach to model selection in ABC (Pudlo et al. 2016; Marin
et al. 2018) by more appropriately targeting the full multivariate
distribution over models in the dependent variable. More details about DRF are
provided later.

Step 3 of Algorithm 1 aims to estimate the posterior distribution of possibly
high-dimensional model parameters at least approximately-well, using a
tractable multivariate copula-based distribution, based on Sklar’s (1959)
theorem.

To explain this theorem in terms of Algorithm 1, let

$\Pi^{\text{ABC}}(\theta\mid
s(x))\equiv\Pi_{x}(\theta_{1},\ldots,\theta_{d})=\Pi_{x}(\theta)$

be the cdf of the approximate posterior distribution (approximate because the
summaries $s$ are not necessarily sufficient), with corresponding univariate
posterior marginal cdfs $\Pi_{x,1}(\theta_{1}),\ldots,\Pi_{x,d}(\theta_{d})$
and pdfs $\pi_{x,1}(\theta_{1}),\ldots,\pi_{x,d}(\theta_{d})$.

Sklar’s (1959) theorem implies that there exists a copula $C$ (cdf) such that
for all $(\theta_{1},\ldots,\theta_{d})\in\mathbb{R}^{d}$, the joint posterior
cdf can be represented by

$\Pi_{x}(\theta_{1},\ldots,\theta_{d})=C\\{\Pi_{x,1}(\theta_{1}),\ldots,\Pi_{x,d}(\theta_{d})\\}.$

If the $\Pi_{x,k}$ (cdfs) are continuous for all $k=1,\ldots,d$ (with
corresponding pdfs $\pi_{x,k}$), then $C$ is unique, with uniform distributed
$\Pi_{x,k}(\theta_{k})\sim$ $\mathcal{U}_{(0,1)}$ over the random $\theta_{k}$
(for $k=1,\ldots,d$) by the probability integral transform, and the copula $C$
is ‘margin-free’ in the sense that it is invariant under increasing
transformations of the margins (Nelsen 2006, Theorem 2.4.3).

Then for an arbitrary continuous multivariate distribution (Okhrin et al.
2017, §13.2), its copula $C$ (cdf) can be determined from the transformation

$C(u_{1},\ldots,u_{d})=\Pi_{x}(\Pi_{x,1}^{-1}(u_{1}),\ldots,\Pi_{x,d}^{-1}(u_{d})),\text{
for }u_{1},\ldots,u_{d}\in[0,1],$

with corresponding copula density (pdf)
$c(u_{1},\ldots,u_{d})=\tfrac{\partial^{d}C(u_{1},\ldots,u_{d})}{\partial
u_{1},\ldots,\partial u_{d}}$ (for $u_{1},\ldots,u_{d}\in[0,1]$), where the
$\Pi_{x,k}^{-1}$ (for $k=1,\ldots,d$) are inverse marginal distribution
functions. Also, the posterior distribution pdf of $\theta$ is given by

$\pi_{x}(\theta_{1},\ldots,\theta_{d})=c\\{\Pi_{x,1}(\theta_{1}),\ldots,\Pi_{x,d}(\theta_{d})\\}{\textstyle\prod\nolimits_{k=1}^{d}}\pi_{x,k}(\theta_{k})\text{
for }\theta_{1},\ldots,\theta_{d}\in\mathbb{R}\cup\\{\infty,-\infty\\}.$

Algorithm 1. Estimate the posterior probability density function (pdf):
---
$\pi_{\widehat{m}}^{\text{ABC}}(\theta\mid
s(x))\propto\pi_{\widehat{m}}(\theta){\textstyle\int}f_{\widehat{m},s}(s(x)+\varepsilon\nu\mid\theta)\mathcal{K}(\nu)\mathrm{d}\nu$,
for the selected (best) model $\widehat{m}$ among $M\geq 1$ models considered.
Inputs: Dataset $x$; number of iterations $N_{M}$; and $M\geq 1$ models
$\\{m\\}_{m=1}^{M}$ for
the data $x$; each model $m$ with exact likelihood $f_{m}(\cdot\mid\theta)$,
prior probability $\pi(m)$,
prior distribution $\pi_{m}(\theta)$ for model parameters
$\theta=(\theta_{k})_{k=1}^{d(m)}\in\Theta^{(m)}\subseteq\mathbb{R}^{d(m)}$,
and candidate summary statistics $s(\cdot)$.
(Step 1) Construct an initial Reference Table.
for $j=1,\ldots,N_{M}$ do (in parallel computing cores if necessary)
| (a) | If $M>1$, draw model index, $m(j)\sim\pi(m)$. If $M=1$, set $m(j)\equiv\widehat{m}\equiv 1$.
| (b) | Draw a parameter sample, $\theta^{(j)}\sim\pi_{m(j)}(\theta)$.
| (c) | Draw a sample dataset, $y^{(j)}\sim f_{m(j)}(x\mid\theta^{(j)})$.
| (d) | Calculate summary statistics, $s(y^{(j)})$.
| (e) | For each discrete $\theta_{k}^{\ast(j)}\in\theta^{(j)}$,
|  | draw $U_{k}^{(j)}\sim$ $\mathcal{U}_{(0,1)}$, then set $\theta_{k}^{(m(j))}\equiv$ $\theta_{k}^{\ast(m(j))}-U_{k}^{(j)}$.
end for
Output Reference Table: $\\{m(j),\theta^{(j)},s(y^{(j)})\\}_{j=1}^{N_{M}}$;
or $\\{\theta^{(j)},s(y^{(j)})\\}_{j=1}^{N}$ if $M=1$, with $N_{M}\equiv N$.
(Step 2) If $M\ >\ 1$, train DRF on $\\{m(j),s(y^{(j)})\\}_{j=1}^{N_{M}}$,
regressing $m$ on $s$
(selected automatically by DRF). Then based on trained DRF, estimate the
posterior probabilities $\widehat{\Pi}_{x}(m)$ of models $\\{m\\}_{m=1}^{M}$
conditional on $s(x)$.
Select model $\widehat{m}\in\\{m\\}_{m=1}^{M}$ with highest posterior
probability.
Reduce Reference Table to $\widehat{m}$: $\
\\{\theta^{(j)},s(y^{(j)})\\}_{j=1}^{N}\equiv$
$\\{\theta^{(j)},s(y^{(j)}):m(j)=\widehat{m}\\}_{j=1}^{N}$.
(Step 3)  Estimate parameters $\\{\nu,\rho,\\{\pi_{k},\Pi_{k}\\}_{k=1}^{d}\\}$
of meta-$t$ posterior pdf:
$\pi_{\widehat{m},\nu,\rho,\pi}^{\text{ABC}}(\theta\mid
x)=c_{\nu,\rho}(u){\textstyle\prod\limits_{k=1}^{d}}\pi_{x,k}(\theta_{k})$
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
=\dfrac{t_{d,\nu,\rho}(T_{\nu}^{-1}(u_{1}),\ldots,T_{\nu}^{-1}(u_{d}))}{{\textstyle\prod\nolimits_{k=1}^{d}}t_{\nu}(T_{\nu}^{-1}(u_{k}))}{\textstyle\prod\limits_{k=1}^{d}}\pi_{x,k}(\theta_{k})$,
for $u\in(0,1)^{d}$;
$c_{\nu,\rho}(u)$ is $t$-copula pdf; $\ \ t_{\nu}$ is pdf ($T_{\nu}$ cdf) of
standard univariate $t$ distribution;
$t_{d,\nu,\rho}$ is the $d$-variate $t$ pdf; with degrees of freedom $\nu$ and
$d\times d$ scale matrix $\rho$;
univariate marginal posterior pdfs $\pi_{x,k}(\theta_{k})$ and cdfs
$u_{k}=\Pi_{x,k}(\theta_{k})$ for $k=1,\ldots,d$.
Find estimates
$\\{\widehat{\nu},\widehat{\rho},\\{\widehat{\pi}_{k},\widehat{\Pi}_{k}\\}_{k=1}^{d}\\}$
using 2-stage semiparametric estimation:
(1) For $k=1,\ldots,d$, train DRF on
$\\{\theta_{k}^{(j)},s(y^{(j)})\\}_{j=1}^{N}$, regressing $\theta_{k}$ on
$s(y)$,
with summary statistics (covariate) selection; use trained DRF to predict
estimates
of posterior cdfs
$\widehat{u}_{k}^{(j)}=\widehat{\Pi}_{x,k}(\theta_{k}^{(j)})$ & pdfs
$\widehat{\pi}_{x,k}(\theta_{k}^{(j)})$ conditional on $s(x)$.
(2) Find: $(\widehat{\nu},\widehat{\rho})=$
$\underset{(\nu,\rho)\in(0,\infty)\times\\{\rho\\}}{\arg\max}{\textstyle\prod\limits_{j=1}^{N}}c_{\nu,\rho}(\widehat{\Pi}_{x,1}(\theta_{1}^{(j)}),\ldots,\widehat{\Pi}_{x,d}(\theta_{d}^{(j)})){\textstyle\prod\limits_{k=1}^{d}}\widehat{\pi}_{x,k}(\theta_{k}^{(j)})$.
(Step 4)  From estimated posterior distribution
$\pi_{\widehat{m},\widehat{\nu},\widehat{\rho},\widehat{\pi}}^{\text{ABC}}(\theta\mid
x)$, obtain posterior:
mean, variance, quantiles of $\theta_{k}$ from DRF cdf estimates
$\widehat{\Pi}_{x,k}$ for $k=1,\ldots,d$;
posterior correlation matrix estimate $\widehat{\rho}$ of $\theta$ from Step
3; and Mode and MLE by:
$\widehat{\text{Mode}}_{x}(\theta)=$
$\underset{j=1,\ldots,N}{\arg\max}\pi_{\widehat{m},\widehat{\nu},\widehat{\rho},\widehat{\pi}}^{\text{ABC}}(\theta^{(j)})$,
MLE $\widehat{\theta}=$
$\underset{j=1,\ldots,N}{\arg\max}\pi_{\widehat{m},\widehat{\nu},\widehat{\rho},\widehat{\pi}}^{\text{ABC}}(\theta^{(j)})/\pi(\theta^{(j)})$;
or by:
$\widehat{\text{Mode}}_{x}(\theta)=\underset{j=1,\ldots,N_{+}}{\arg\max}\pi_{\widehat{m},\widehat{\nu},\widehat{\rho},\widehat{\pi}}^{\text{ABC}}(\theta_{+}^{(j)})$,
$\widehat{\theta}=\underset{j=1,\ldots,N_{+}}{\arg\max}\pi_{\widehat{m},\widehat{\nu},\widehat{\rho},\widehat{\pi}}^{\text{ABC}}(\theta_{+}^{(j)})/\pi(\theta_{+}^{(j)})$,
given $N_{+}\ $additional draws $\\{\theta_{+}^{(j)}\\}_{j=1}^{N_{+}}$
$\overset{\text{iid}}{\sim}$
$\pi_{\widehat{m},\widehat{\nu},\widehat{\rho},\widehat{\pi}}^{\text{ABC}}$.

As mentioned, the ABC framework (Algorithm 1) employs a novel multivariate
meta-$t$ (Fang et al. 2002) copula-based distribution approximation of the
posterior, which allows for multivariate skewness by using skewed marginals
which are modeled flexibly and accurately estimated by DRFs, respectively.
This multivariate copula-based distribution has marginals which are
(covariate-dependent) regressions and a multivariate $t$ parametric copula
that is independent of covariates, extending previous related copula models
(Pitt et al. 2006; Song, 2000) through the use of DRF. The meta-$t$
distribution is defined by a $t$ copula cdf $C_{\nu,\rho}$ (copula density
pdf, $c_{\nu,\rho}$) with degrees of freedom $\nu$ and $d\times d$ correlation
matrix parameters $\rho$.

As an aside, when the random vector $\theta$ contains any discrete variables,
the joint cdf $C$ is not identifiable under Sklar’s theorem (e.g., Geenens,
2020). Therefore, to maintain direct use of Sklar’s theorem using the
$t$-copula, Step 1(e) of Algorithm 1 can routinely apply jittering (Denuit et
al. 2005; Madsen & Fang, 2011) to continuize each (of any) discrete model
parameter (integer-valued without loss of generality), while ensuring that no
information is lost and preserving the same dependence relationship among all
model parameters. Each (of any) discrete (integer) model parameter
$\theta_{k}^{\ast}\in\theta$ with posterior cdf $\Pi_{x,k}^{\ast}(\theta_{k})$
is continuized (“jittered”) to continuous parameter $\theta_{k}$ by
$\theta_{k}=\theta_{k}^{\ast}-U_{k}$ using uniform random variable $U_{k}\sim$
$\mathcal{U}_{(0,1)}$ (with $U_{k}\perp\theta_{k}^{\ast}$ and $U_{k}\perp
U_{l}$ for $k\neq l$), so that $\theta_{k}$ has posterior cdf

$\Pi_{x,k}(\theta_{k})=\Pi_{x,k}^{\ast}(\lfloor\theta_{k}\rfloor)+(\theta_{k}-\lfloor\theta_{k}\rfloor)\Pr(\vartheta_{k}^{\ast}=\lfloor\theta_{k}+1\rfloor)$

and pdf
$\pi_{x,k}(\theta_{k})=\Pr(\vartheta_{k}^{\ast}=\lfloor\theta_{k}+1\rfloor)$,
with $\lfloor z\rfloor=$ floor$(z)$. The parameters of $\Pi_{x,k}$ and
$\pi_{x,k}$ are exactly those of $\Pi_{x,k}^{\ast}$, and $\theta_{k}^{\ast}$
can be recovered from $\theta_{k}$ as
$\theta_{k}^{\ast}=\lfloor\theta_{k}+1\rfloor$.

Step 3 of Algorithm 1 performs two-stage semiparametric estimation (Genest et
al. 1995) of the meta-$t$ posterior distribution parameters,
$\\{\nu,\rho,\\{\pi_{k},\Pi_{k}\\}_{k=1}^{d}\\}$. The first stage calculates
values of the nonparametric marginal posterior cdfs and pdfs,
$\\{\widehat{u}_{k}^{(j)}=\widehat{\Pi}_{x,k}(\theta_{k}^{(j)}),\widehat{\pi}_{x,k}(\theta_{k}^{(j)})\\}_{j=1}^{N},$
conditionally on $s(x)$, from DRFs trained on reference subtables
$\\{\theta_{k}^{(j)},s(y^{(j)})\\}_{j=1}^{N}$ for $k=1,\ldots,d$ (resp.),
while performing automatic summary statistics (covariate) selection. The
second stage implements an Expectation-Maximization (EM) algorithm (Hintz et
al 2022a,b) to calculate the MLE $(\widehat{\nu},\widehat{\rho})$ of the
parameters of the meta-$t$ pdf based on these estimated cdfs and pdfs.

More specifically, DRF is trained (estimated) on the reference table
$\\{\theta_{k}^{(j)},s(y^{(j)})\\}_{j=1}^{N}$ (the training dataset), such
that prior parameter samples $\theta_{k}^{(j)}$ (dependent variable) are
regressed on the summary statistics (e.g. MPLEs) covariate vectors
$s(y^{(j)}),$ while targeting the conditional distribution of the dependent
variable, for each model parameters indexed by $k=1,\ldots,d$. Then based on
each of these $d$ trained DRFs, the posterior distribution (cdfs $\Pi_{x,k}$
and pdfs $\pi_{x,k}$, for $k=1,\ldots,d$) can be accurately predicted
(estimated) (Ćevid et al. 2022), conditionally on the summary statistics of
the original observed dataset, $s(x)$; all while performing automatic summary
statistics (variable) selection from a potentially larger set of summary
statistics, and accounting for the uncertainty in variable selection, without
requiring the user to pre-select the subset of summary statistics relevant to
the intractable-likelihood model under consideration. DRF is a weighted
nearest neighbor method (Lin & Jeon, 2006) which performs locally-adaptive
estimation of the conditional distribution through the aggregation of
dependent variable predictions of an ensemble of randomized flexible
Classification of Regression Trees (CARTs) (Breiman et al. 1984),
respectively, such that each of these CARTs is estimated from a random
subsample of training dataset (drawn without replacement), and that each level
of the CART’s binary tree is constructed by splitting the training data points
$\\{\theta_{k}^{(j)},s(y^{(j)})\\}_{j=1}^{N}$ based on a covariate
$s_{k}(y^{(j)})\in s(y^{(j)})$ and its split point $c$ both chosen in such a
way that the distribution of the dependent responses $\theta_{k}^{(j)}$ for
which $s_{k}(y^{(j)})\leq c$ differs the most compared to the distribution for
which $s_{k}(y^{(j)})>c$, according to the Maximal Mean Discrepancy (MMD)
(Gretton et al. 2007), a quickly-computable two-sample test statistic that can
detect a wide variety of distributional changes. This way, the dependent
variable distribution in each of the resulting leaf nodes is as homogeneous as
possible, to define neighborhoods of relevant training data points for every
covariate value $s$. Repeating this many times with randomization induces a
weighting function $w_{s(x)}(s(y^{(j)}))$ quantifying the relevance of each
training covariate datapoint $s(y^{(j)})$ (for $j=1,\ldots,N$) for a given
test point $s(x)$. The DRF estimates of the posterior conditional cdfs are
given by:

$\displaystyle\widehat{\Pi}_{x,k}(\theta_{k})$
$\displaystyle={\displaystyle\sum\limits_{j=1}^{N}}w_{s(x)}\\{s(y^{(j)})\\}\mathbf{1}(\theta_{k}^{(j)}\leq\theta_{k})$
(2.1a)
$\displaystyle={\displaystyle\sum\limits_{j=1}^{N}}\left\\{\dfrac{1}{B}{\displaystyle\sum\limits_{b=1}^{B}}\frac{\mathbf{1}[s(y^{(j)})\in\mathcal{L}_{b}\\{s(x)\\}]}{\left|\mathcal{L}_{b}\\{s(x)\\}\right|}\right\\}\mathbf{1}(\theta_{k}^{(j)}\leq\theta_{k}),\text{
}$ (2.1b) for $k=1,\ldots,d$, obtained from $N$ (training) points of the
Reference table $\\{\theta_{k}^{(j)},s(y^{(j)})\\}_{j=1}^{N}$; where
$w_{s(x)}(s(y^{(j)}))$ are positive weights with
${\textstyle\sum\nolimits_{j=1}^{N}}w_{s(x)}\\{s(y^{(j)})\\}=1$; and $B$ is
the number of randomized CARTs ($T_{b}$ for $b=1,\ldots,B$) of the ensemble;
$\mathcal{L}_{b}\\{s(x)\\}$ and $\left|\mathcal{L}_{b}\\{s(x)\\}\right|$ are
respectively the set and number of the training data points
$(\theta^{(j)},s(y^{(j)}))$ which end up in the same leaf as $s(x)$ in CART
tree $T_{b}$; and $\mathbf{1}(\cdot)$ is the indicator function. The
corresponding DRF pdf estimates (for $k=1,\ldots,d$) can be found by the
following empirical histogram density estimator:
$\widehat{\pi}_{x,k}(\theta_{k})={\displaystyle\sum\limits_{l=1}^{N}}\frac{w_{s(x)}\\{s(y_{(l)})\\}}{\theta_{k,(l)}-\theta_{k,(l-1)}}\mathbf{1}(\theta_{k}\in(\theta_{k,(l-1)},\theta_{k,(l)}]),$
(2.2)

where $\theta_{k,(1)},\ldots,\theta_{k,(N)}$ (with
$\theta_{k,(0)}\equiv\theta_{k,(1)}$) are the order statistics of
$\theta_{k}^{(1)},\ldots,\theta_{k}^{(N)}$, these order statistics
corresponding to weights
$w_{s(x)}\\{s(y_{(1)})\\},\ldots,w_{s(x)}\\{s(y_{(N)})\\}$.

Extensive simulation studies (Ćevid et al. 2022) showed that DRF performs well
and outperforms other machine learning methods in terms of prediction
accuracy, for a wide range of sample size and problem dimensionality,
especially in problems where $p$ (number of covariates) is large and $d$
(dimensionality of the dependent variable) is small to moderately large,
without the need for further tuning or involved numerical optimization. DRF
training and predictions can compute fast using the drf R package (Michel &
Cevid, 2021) with minimal tuning parameters. By default, a DRF is estimated
from an ensemble of $B=2,000$ CARTs, with every tree constructed (randomized)
from a random subset being 50% of the size of the training set, with a target
of 5 for the minimum number of observations in each CART tree leaf (Ćevid et
al. 2022, p.33). The estimated (induced) weighting function of DRF, for each
of the given model parameters $\theta_{k}^{(j)}$ (for $k=1,\ldots,d$), is used
to estimate a parametric meta-$t$ multivariate copula-based distribution for
these dependent variable observations, after having nonparametrically adjusted
for the covariates $s(y^{(j)})$ (for $j=1,\ldots,N$) (in the spirit of Bickel
et al. 1993). Optional Step 2 of Algorithm 1, used when a set of $M>1$
Bayesian models is considered, performs model selection based on DRF estimates
of the $M$ posterior model probabilities, conditionally on observed data
summaries $s(x)$.

Step 4 of Algorithm 1, with the $d$ DRFs estimated in Step 3, uses the weights
$w_{s(x)}$ and marginal posterior cdf estimates (2.1) to calculate, for each
model parameter $\theta_{k}$ for $k=1,\ldots,d$, estimates of marginal
posterior expectations, including estimates of the posterior mean
($\widehat{\mathbb{E}}$), variance ($\widehat{\mathbb{V}}$), and quantiles
($\widehat{Q}$):

$\displaystyle\widehat{\mathbb{E}}_{x,k}(\theta_{k})$
$\displaystyle={\textstyle\sum\nolimits_{j=1}^{N}}\theta_{k}^{(j)}w_{s(x)}\\{s(y^{(j)})\\},$
$\displaystyle\widehat{\mathbb{V}}_{x,k}(\theta_{k})$
$\displaystyle={\textstyle\sum\nolimits_{j=1}^{N}}\\{\theta_{k}^{(j)}-\widehat{\mathbb{E}}_{x,k}(\theta_{k})\\}^{2}w_{s(x)}\\{s(y^{(j)})\\},$
$\displaystyle\widehat{Q}_{x,k}(u)$
$\displaystyle=\widehat{\Pi}_{x,k}^{-1}(u),\text{ for }k=1,\ldots,d\text{ and
}u\in(0,1).$

Finally, Step 4 calculates estimates of the multivariate posterior mode,
$\widehat{\text{Mode}}_{x}(\theta)$, and the MLE $\widehat{\theta}$ of the
parameters $\theta$ of the (selected) model. The new multivariate mode
estimator maximizes a (posterior) density estimate (or finds the MLE by
maximizing a ratio of the posterior density estimate to the prior density)
over a sample of points, being the parameter samples $\theta^{(j)}$ from the
Reference Table, or if necessary, an extra set of parameter samples
$\theta_{+}^{(j)}$ drawn from the meta-$t$ posterior density estimate. The
estimation of the posterior mode and MLE is done while performing automatic
summary statistics (variable) selection via DRFs for the meta-$t$ marginals.

##  3 Numerical Examples

In §3.1, we present results of simulation studies for models with analytically
available and directly calculable posterior distributions, allowing direct
evaluation of how well the ABC Algorithm 1 estimates the true posterior
distribution. There, we evaluate results for a typical strategy where the size
($n_{\text{sim}}$) of the dataset $y$ simulated from the model in each
iteration of ABC Algorithm 1 matches the size ($n$) of the dataset $x$ being
analyzed (i.e., $c=n_{\text{sim}}$), and a strategy where the dataset
simulated from the model per algorithm iteration is smaller than the dataset
being analyzed, i.e., $n_{\text{sim}}<n$, while the summary statistics
employed take on values which have meaning that is invariant to the size of
the dataset. The latter strategy may be considered in situations where a large
or massive dataset is being analyzed and/or simulating from the model is
computationally costly, where the approach of simulating datasets of smaller
size $n_{\text{sim}}<n$ help reduce the computational cost. Further related
discussion is provided below.

In §3.2-§3.8, we consider ABC inference for models defined by intractable
likelihoods, from which it is not easy to draw samples from the model
likelihood and/or to calculate point-estimators from, especially when the
dataset sample size is large, as §1 mentioned. Specifically, we consider ERGM
and mechanistic models for network datasets, reviewed in §3.2, to help set the
scenes in §3.3-§3.8. These models are applied in simulation studies in
§3.3-§3.5, and applied to analyze real-life massive network datasets in
§3.6-§3.8. In such large network scenarios, it is prohibitively costly and not
practically feasible to simulate many networks of the same size (number of
nodes) as the original network dataset over iterations in an ABC or synthetic
likelihood algorithm, especially considering that already it can be
prohibitively costly to simulate one network from ERGM or mechanistic network
models (e.g., ERGM given parameters simulates a network using MCMC) for a
network of sufficiently large size (number of nodes), let alone to compute
point-estimates (e.g., MLE) of the parameters of the given network model.
These issues, in the context of ABC Algorithm 1, are addressed by a strategy
that simulates network datasets of a smaller size compared to the size of the
network dataset under statistical analysis, while using network summary
statistics (calculated on the observed dataset and each simulated dataset)
that take on values that have the same meaning and are invariant to the size
of the network dataset(s) being analyzed. These network size invariant summary
statistics are reviewed in §3.2. The simulation studies in §3.3-§3.5 evaluate
the posterior distribution estimation accuracy of ABC Algorithm 1 based on
this strategy, which is also used in the real data applications of §3.6-§3.8.

Each application of ABC Algorithm 1 in §3 is based on a Reference Table of
size N = 10,000 samples, and on using R code and packages (Michel & Ćevid
2021; Hintz et al. 2022a, Handcock et al. 2023; Krivitsky, 2022; Csardi et al.
2024).

###  3.1 Simulation Study: Calculable Posterior Distributions

Here we study the accuracy of ABC Algorithm 1 to estimates the directly-
calculable true posterior distribution of a Poisson$(\theta)$ model and a two-
component scale normal mixture model with common mean parameter $\mu$, which
has been studied in previous research on ABC and synthetic likelihood methods.

The Poisson$(\lambda)$ distribution, is defined by a tractable likelihood
probability mass function (p.m.f.)
$f(X_{1}=x_{1}\mid\lambda)=\lambda^{x_{1}}\exp(-\lambda)/x_{1}!$, is assigned
a gamma prior distribution $\mathcal{G}(\lambda\mid\frac{1}{2},0.1)$ (shape
$\frac{1}{2}$ and rate $0.1$), and exact gamma posterior distribution
$\lambda\mid
x_{1}\sim\mathcal{G}(\lambda\mid\frac{1}{2}+{\textstyle\sum\nolimits_{i=1}^{n}}x_{i,1},0.1+n)$,
with the sample mean
$\overline{x}_{1}=\tfrac{1}{n}{\textstyle\sum\nolimits_{i=1}^{n}}x_{i,1}$ a
sufficient summary statistic for $\lambda$. The scale normal mixture model
$f(X_{2}=x_{2}\mid\mu)=\frac{1}{2}\mathcal{N}(x_{2}\mid\mu,1)+\frac{1}{2}\mathcal{N}(x_{2}\mid\mu,0.01)$,
the normal distribution $\mathcal{N}(x_{2}\mid\mu,z+(1-z)0.01)$ with latent
covariate $z\sim$ Bern$(1/2)$, is assigned a uniform prior distribution
$\mu\sim\mathcal{U}_{(-10,10)}$, yielding the posterior distribution
$\pi(\mu\mid
x_{2})\propto\frac{1}{2}\mathcal{N}(\mu\mid\overline{x}_{2},1)+\frac{1}{2}\mathcal{N}(\mu\mid\overline{x}_{2},0.01)$
(effectively, with equality up to machine precision) with summary statistic
the sample mean
$\overline{x}_{2}=\tfrac{1}{n}{\textstyle\sum\nolimits_{i=1}^{n}}x_{i,2}$. The
Poisson and scale normal mixture distributions are specified to define a
bivariate distribution model, respectively, with parameters assigned the
aforementioned independent prior distributions, and posterior distributions as
marginal distributions of $\theta=(\lambda,\mu)$. This enables ABC Algorithm 1
to estimate a bivariate posterior mode and MLE of $\theta$.

For the simulation study, a total of 10 replications of $n=100$ bivariate
observations $\\{(x_{i,1},x_{i,2})\\}_{i=1}^{n=100}$ were simulated from the
bivariate distribution model with true data-generating parameters
$\theta=(\lambda=3,\mu=0)$, for each of three conditions defined by
$n_{\text{sim}}$ $=100,50,$ and $33$ simulated dataset samples sizes used for
each iteration of ABC Algorithm 1. Intuitively, the sufficient statistic
$\overline{x}_{1}$ for the Poisson distribution parameter $\lambda$ is
sufficient for a dataset of size $n_{\text{sim}}$ $=100$, and only
approximately sufficient for the smaller simulation sample sizes
$n_{\text{sim}}$ $=50$ and $33$. Similarly for the statistic
$\overline{x}_{2}$ summarizing the parameter $\mu$ of the scale normal mixture
model.

Table 1: Results of the simulation study for the Poisson$(\lambda)$ and
NormalMixture$(\mu)$ model.
---
Mean (standard deviation) of posterior estimate over 10 replicas.
| $\lambda$ | $\mu$ | $\lambda$ | $\mu$ | $\lambda$ | $\mu$
n | 100 | 100 | 100 | 100 | 100 | 100
n${}_{\text{sim}}$ | 100 | 100 | 50 | 50 | 33 | 33
true $(\lambda,\mu)$ | 3.00 | 0.00 | 3.00 | 0.00 | 3.00 | 0.00
ABC mean | 3.02 (0.15) | 0.01 (0.10) | 3.09 (0.17) | -0.02 (0.10) | 3.03 (0.13) | 0.01 (0.07)
Exact mean | 3.02 (0.15) | 0.01 (0.10) | 3.08 (0.17) | -0.02 (0.09) | 2.99 (0.13) | 0.02 (0.08)
ABC med. | 3.02 (0.15) | 0.01 (0.10) | 3.09 (0.16) | -0.02 (0.09) | 3.03 (0.13) | 0.01 (0.07)
Exact med. | 3.02 (0.15) | 0.01 (0.10) | 3.07 (0.17) | -0.02 (0.09) | 2.98 (0.13) | 0.02 (0.08)
ABC mode | 3.06 (0.24) | -0.14 (0.44) | 3.07 (0.27) | -0.16 (0.27) | 3.09 (0.27) | -0.12 (0.26)
Exact mode | 3.01 (0.15) | 0.01 (0.10) | 3.07 (0.17) | -0.02 (0.09) | 2.98 (0.13) | 0.02 (0.08)
ABC MLE | 3.06 (0.24) | -0.14 (0.44) | 3.07 (0.27) | -0.16 (0.27) | 3.09 (0.27) | -0.12 (0.26)
Exact MLE | 3.02 (0.15) | 0.01 (0.10) | 3.08 (0.17) | -0.02 (0.09) | 2.99 (0.13) | 0.02 (0.08)
ABC s.d. | 0.20 (0.01) | 0.22 (0.07) | 0.26 (0.02) | 0.24 (0.08) | 0.32 (0.03) | 0.25 (0.11)
Exact s.d. | 0.17 (0.00) | 0.71 (0.00) | 0.18 (0.00) | 0.71 (0.00) | 0.17 (0.00) | 0.71 (0.00)
KS distance | 0.09 (0.02) | 0.20 (0.01) | 0.14 (0.03) | 0.18 (0.01) | 0.20 (0.04) | 0.17 (0.01)
KS test stat* | 0.88 (2.64) | 2.18 (6.58) | 1.74 (5.31) | 2.09 (6.35) | 2.64 (8.10) | 1.79 (5.38)
Copula | d.f. | Correlation | d.f. | Correlation | d.f. | Correlation
d.f. & corr. | 7.89 (7.27) | 0.19 (0.42) | 10.74 (8.60) | 0.02 (0.40) | 21.28 (11.56) | 0.03 (0.16)
* Kolmogorov-Smirnov (KS) test statistic 95% critical value equals 1.358.
MAE (MSE) of posterior estimate over 10 replicas.
| $\theta$ | $\mu$ | $\theta$ | $\mu$ | $\theta$ | $\mu$
n | 100 | 100 | 100 | 100 | 100 | 100
n${}_{\text{sim}}$ | 100 | 100 | 50 | 50 | 33 | 33
ABC mean | 0.10 (0.02) | 0.09 (0.01) | 0.14 (0.04) | 0.07 (0.01) | 0.09 (0.02) | 0.06 (0.00)
Exact mean | 0.11 (0.02) | 0.08 (0.01) | 0.15 (0.03) | 0.06 (0.01) | 0.09 (0.02) | 0.07 (0.01)
ABC med. | 0.11 (0.02) | 0.08 (0.01) | 0.14 (0.03) | 0.07 (0.01) | 0.10 (0.02) | 0.06 (0.00)
Exact med. | 0.11 (0.02) | 0.08 (0.01) | 0.15 (0.03) | 0.06 (0.01) | 0.09 (0.02) | 0.07 (0.01)
ABC mode | 0.20 (0.06) | 0.30 (0.19) | 0.17 (0.07) | 0.20 (0.09) | 0.21 (0.07) | 0.16 (0.08)
Exact mode | 0.11 (0.02) | 0.08 (0.01) | 0.15 (0.03) | 0.06 (0.01) | 0.09 (0.02) | 0.07 (0.01)
ABC MLE | 0.20 (0.06) | 0.30 (0.19) | 0.17 (0.07) | 0.20 (0.09) | 0.21 (0.07) | 0.16 (0.08)
Exact MLE | 0.11 (0.02) | 0.08 (0.01) | 0.15 (0.03) | 0.06 (0.01) | 0.09 (0.02) | 0.07 (0.01)

For each of the two parameters of the bivariate model and each of the three
simulation conditions ($n_{\text{sim}}$ $=100,50,$ and $33$), the accuracy of
the marginal posterior distribution estimate of ABC Algorithm 1, relative to
the known exact marginal posterior distribution of the parameter (see above)
is measured by the weighted one-sample (two-tailed) Kolmogorov-Smirnov (KS)
distance and corresponding test statistic (Monohan 2022, p. 358) based on the
estimated DRF weights, with the test statistic known to have 95th percentile
critical value of 1.358. Also, for ABC Algorithm 1, its posterior mean,
median, and mode estimate of each model parameter is measured by median
absolute error (MAE) and by mean squared error (MSE) which are, respectively,
the absolute distance (for MAE) and the squared distance (for MSE) between the
given estimate and the true data-generating model parameter, averaged over the
10 simulated data replicas. Table 1 presents the results of the simulation
study, which shows that ABC Algorithm 1 produced posterior distribution and
point estimates comparable to the respective exact values for all conditions
$n_{\text{sim}}$ $=100,$ $50$ and $33$, with the results for the point
estimates (according to MAE and MSE) being more favorable than the results for
the full posterior distribution (according to the KS statistics). This shows
that simulating datasets from the given model of size $n_{\text{sim}}$ smaller
than the size $n$ of the observed dataset being analyzed can be a reasonable
approach for obtaining accurate point estimates using ABC Algorithm 1.

###  3.2 Overview of Network Data Modeling

Network science which studies systems and phenomena governing networks,
representing relational data which shape the connections of the world and are
ubiquitous and focal discussion points in everyday life and the sciences. A
network is comprised of a set of nodes (i.e., vertices, actors, individuals,
or sites, etc.) which can have connections (i.e., edges, ties, lines, dyadic
relationships, links, or bonds, etc.) with other nodes in the network. For
example, a social network may consist of friendship ties among persons
(nodes). A directed network treats each pair of nodes $(i,j)$ and $(j,i)$ as
distinct, whereas an undirected network does not. Each possible edge tie
connecting any given node pair in a network may be (unweighted) binary-valued,
or (weighted) count-, ordinal-, categorical-, text-, ranking-, continuous
valued; and/or vector valued; possibly as part of a multi-layered (or
multiplex, multi-item, or multi-relational) network representing multiple
relationships among the same set of nodes. Network summary statistics include
network size (number of nodes, $n$), edge count and density (ratio of edge
count to maximum edge count in network), node degree (a node’s outdegree is
the number of edges starting from it, indegree is the number of edges going
into it, while in an undirected network, the outdegree and indegree coincide
and equal to degree), the degree distribution over nodes and its geometric
weighted form (Snijders et al. 2006), degree assortativity (Newman, 2003),
counts of 2-stars and triangles, and global (Wasserman & Faust, 1994) and
average local (Watts & Strogatz, 1998) clustering coefficients. The presence
of a dyadic network tie can depend on the presence of other dyadic ties, nodal
attributes, the network’s evolving dynamic structure and dyadic interactions
over time, or on other factors; and the network’s structure itself may
influence or predict individuals’ attributes and behaviors.

Network data are collected and analyzed in various scientific fields,
including the social sciences (e.g., social or friendship networks, marriage,
sexual partnerships); academic paper networks (collaboration or citation
networks); neuroscience (human brain networks and interactions); political
science (international relations, wars between countries, insurgencies,
terrorist networks, strategic alliances and friendships); education
(multilevel student relation data, item response data); economics (financial
markets, economic trading or resource distribution networks); epidemiology
(disease spread dynamics, HIV, COVID-19); physics (Ising models, finding
parsimonious mechanisms of network formation); biology and other natural
sciences (protein-protein, molecular interaction, metabolic, cellular,
genetic, ecological, and food web networks); artificial intelligence and
machine learning (finding missing link in a business or a terrorist network,
recommender systems, Netflix, neural networks); spatial statistics (discrete
Markov random fields); traffic systems and transportation networks (roads,
railways, airplanes); cooperation in an organization (advice giving in an
office, identity disambiguation, business organization analysis);
communication patterns and networks (detecting latent terrorist cells);
telecommunications (mobile phone networks); and computer science and networks
(e-mail, internet, blogs, web, Facebook, Twitter(X), LinkedIn, information
spread, viral marketing, gossiping, dating networks, blockchain network, virus
propagation). Network science has seen many reviews due to its longtime
scientific and societal impacts (Wasserman & Faust, 1994; Newman 2003;
Goldenberg et al. 2010; Snijders, 2011; Pastor-Satorras & Vespignani, 2004;
Raval & Ray, 2013; Newman 2018; Horvat & Zweig, 2018; Brodka & Kazienko 2018;
Ghafouri & Khasteh, 2020; Loyal & Chen 2020; Hammoud & Kramer 2020; Kinsley et
al. 2020; Caimo & Gollini, 2022). Meanwhile, the high dimensionality of
network data poses challenges to modern statistical computing methods for
network data analysis modeling.

Statistical and mechanistic network models are two prominent paradigms for
analyzing network data. A statistical network model specifies an explicit
likelihood function for the given network dataset, available in closed-form up
to a normalizing constant, making standard statistical inference tools
generally available for these models, e.g., for parameter estimation and model
selection. Statistical network models include the popular large family of
Exponential Random Graph models (ERGMs)
(Frank & Strauss 1986; Wasserman & Pattison, 1996; Lusher et al. 2013; Harris
2013; Schweinberger et al. 2020; Caimo & Gollini, 2022) which uses observable
network configurations (e.g., k-stars, triangles) as sufficient statistics. An
ERGM is defined by a likelihood with normalizing constant which is intractable
for a network with more than a handful of nodes $n$. The ERGM for a binary
undirected or directed network is as follows:

$f(X=x\mid\theta)=\dfrac{\exp\left\\{\eta(\theta)^{\intercal}h(x;z)\right\\}}{{\textstyle\sum\nolimits_{x^{\prime}\in\mathcal{X}_{n}}}\exp\left\\{\eta(\theta)^{\intercal}h(x^{\prime};z)\right\\}},$
(3.1)

with parameter vector $\theta\in\Theta\subseteq\mathbb{R}^{q}$ and its mapping
to canonical parameters $\eta:\Theta\rightarrow\mathbb{R}^{d}$ (with (3.1) a
linear ERGM if $\eta(\theta)\equiv\theta$, and is otherwise a curved ERGM);
sufficient statistics $h\in\mathbb{R}^{d}$ which are possibly dependent on any
covariates $z$; and for undirected binary networks $X$, the space
$\mathcal{X}_{n}$ of allowable networks on $n$ nodes is
$\mathcal{X}_{n}=\\{x\in\mathbb{R}^{n\times
n},x_{i,j}=x_{j,i}\in\\{0,1\\},x_{i,i}=0\\}$, and for directed binary networks
$X$, the space $\mathcal{X}_{n}$ is the same but without the restrictions
$x_{i,j}=x_{j,i}$. The ERGM likelihood (3.1) is intractable for a binary
network with more than a handful of nodes $n$, as the normalizing constant of
the model’s likelihood is a sum of $2^{n(n-1)/2}$ terms (or $2^{n(n-1)}$ terms
for a directed network) over the sample space of allowable networks. The basic
ERGM form (3.1) allows it to be straightforwardly extended to a general ERGM
representation for multi-layered weighted (valued) networks shown later in
§3.8, which combines certain ERGMs (Krivitsky 2012; Krivitsky et al. 2020;
Caimo & Gollini, 2022) and encompasses ERGMs: for binary networks; with
size($n$) offset adjusted (invariant) parameterizations (Krivitsky et al 2011;
Stewart et al. 2019); with nodal random effects (Thiemichen et al. 2016); and
for valued (weighted) (Krivitsky 2012; Caimo & Gollini, 2022), multilayered
(Krivitsky et al. 2020), and dynamic networks (Hanneke et al. 2010; Krivitsky
& Handcock, 2014; Lee et al. 2020).

A mechanistic network model is able to incorporate domain scientific knowledge
into generative mechanisms for network formation, enabling researchers to
investigate complex systems using both simulation and analytical techniques. A
mechanistic network model is an algorithmic description of network formation
defined by a small number of domain-specific stochastic microscopic
mechanistic rules by which a network grows and evolves over time, which are
hypothesized and informed by understanding and related critical elements of
the given scientific problem. A typical mechanistic model generates a network
by starting with a small seed network, and then grows the network one node at
a time according to the model’s generative mechanism until some stopping
condition is met, e.g., until a requisite number of nodes is reached, which
may possibly number in the millions. There are easily hundreds of mechanistic
network models, which originated in physics. For a long time, they were
essentially the only type of network models formulated and studied using both
mathematical and computer simulation methods. Mechanistic network models
include the model of (Price 1965, 1976) used to study citation patterns of
scientific papers, and the Barabasi and Albert (1999) model, which introduce
preferential attachment rules for directed and undirected networks,
respectively; the Watts & Strogatz (1998) model, which produces random graphs
with small-world properties, including short average path lengths and high
clustering; the duplication-divergence class of models which describe the
evolution of protein-protein interaction networks (where each node is a
protein in an organism and two proteins are connected if they are involved in
a chemical reaction together), including the Duplication-Mutation-
Complementation (DMC) model (Vazquez et al. 2003; Vazquez 2003) and the
Duplication-Mutation-Random (DMR) model (Sole et al. 2002, Pastor-Satorras et
al. 2003); the KM model (Kretzschmar & Morris 1996; Morris & Kretzschmar,
1997) for modeling HIV epidemics (Morris et al. 2007; Palombi et al. 2012);
and other mechanistic models (Klemm & Eguiluz, 2002; Kumpula et al. 2007;
Procopio et al. 2023). This paper studies a few notable mechanistic network
models, described further in §3.4-§3.5.

We will consider applications of ABC Algorithm 1 to ERGM and some mechanistic
network models through simulation studies §3.3-§3.5 and analyses and real-life
large networks in §3.6-§3.8. Algorithm 1 can be used to analyze a network
dataset $x$ using $M\geq 1$ ERGM and/or mechanistic network models, based on
$K$ scalar and/or vector network statistics, $g_{k}(x)$, for $k=1,\ldots,K$.
It was mentioned that for the analysis of a large network, it is
computationally prohibitive to simulate a very large network from the ERGM or
mechanistic network model, and that this issue can be addressed by a strategy
which simulates network datasets of a smaller size compared to the size of the
network dataset under statistical analysis, while using network summary
statistics (calculated on the observed dataset and each simulated dataset)
that take on values that are invariant to the size (number of nodes) of the
network dataset being analyzed.

In terms of ABC Algorithm 1, we can specify computationally efficient network
summary statistics $s$ which can be used to directly compare between two
networks that may have different node sets, size (number of nodes), and
densities. In other words, summaries which allow for Unknown Node-
Correspondence (UNC) network comparisons (Tantardini et al. 2019, p.2).
Specifying summary statistics vector $s$ by computationally efficient UNC
summary statistics allows Step 1(c) to simulate and compute the MPLEs of a
network of a smaller size $n_{\text{sim}}$, and to directly compare these
summaries with the corresponding summaries of a larger observed network
dataset $x$ (of size $n$) being analyzed, for each of the $N_{M}$ sampling
iterations of ABC Algorithm 1. This lowers computational cost, compared to
comparing doing such summary comparisons based on simulating networks of the
same size ($n$) as the large observed network in each sampling iteration.
Computationally efficient UNC summary statistics include global and average
local clustering coefficients (based on triangle counts) (Pržulj et al. 2004;
Yaverouğlu et al. 2014; Yaverouğlu et al. 2015; Faisal et al. 2017), degree
distribution including its mean and variance or geometric form, degree
assortativity (propensity for similar nodes to connect), and diameter (average
of shortest path lengths over all pairs of nodes in a network) (Yaverouğlu et
al. 2014; Faisal et al. 2017).

A novel contribution is that for either the ERGM or a mechanistic network
model, the UNC network summary statistics vector $s(\cdot)$ may also include
Maximum Pseudo Likelihood Estimates (MPLEs) of $K$ ERGMs defined by sufficient
statistics $g_{k}(x)$, respectively for $k=1,\ldots,K$, with these MPLEs
adjusted by network size by specifying in each ERGM an offset term $|x|$ (edge
count of network $x$) with fixed coefficient $\log\genfrac{(}{)}{}{1}{1}{n}$
(Krivitsky et al. 2011). The MPLE (e.g., Schmid & Hunter, 2023) was introduced
for lattice models (Besag, 1974), and developed to estimate ERGM parameters
(Strauss & Ikeda, 1990; Frank & Strauss, 1986) because the exact MLE of the
ERGM is intractable for a binary network with more than a handful of nodes
$n$, as mentioned above. Specifically, for example, if the observed network
dataset $x$ is an undirected binary network, $x=(x_{i,j}\in\\{0,1\\})_{n\times
n}$ on $n$ nodes, described by scalar or vector valued statistics $g_{k}(x)$
for $k=1,\ldots,K$ (where each $g_{k}$ may depend on covariates $z_{k}$), then
the summary statistics $s$ of $x$ (or $y$) can be specified as
$s(x)=(\widehat{\beta}_{\text{MPLE},1}^{(x)},\ldots,\widehat{\beta}_{\text{MPLE},K}^{(x)})^{\intercal}$,
where $\widehat{\beta}_{\text{MPLE},k}^{(x)}$ is the MPLE of the ERGM (3.1)
using sufficient statistics $g_{k}(x)$ and offset term the network edge count
$|x|$ with fixed coefficient $\log(1/n)$; with each MPLE summary statistic
$\widehat{\beta}_{\text{MPLE},k}$ (for $k=1,\ldots,K$) obtained by maximizing
a logistic regression likelihood, as follows:

$\widehat{\beta}_{\text{MPLE},k}^{(x)}=\underset{\beta_{k}\in\mathbb{R}^{\underline{q}_{k}}}{\arg\max}\underset{1\leq
i<j\leq
n}{{\textstyle\sum}{\textstyle\sum}}\left(\begin{array}[c]{l}x_{i,j}\\{(\eta_{k}(\beta_{k})^{\intercal},\log\genfrac{(}{)}{}{1}{1}{n})^{\intercal}\Delta_{i,j,k}\\}\\\
\lx@intercol\hfil-\log[1+\exp\\{(\eta_{k}(\beta_{k})^{\intercal},\log\genfrac{(}{)}{}{1}{1}{n})^{\intercal}\Delta_{i,j,k}\\}]\hfil\lx@intercol\end{array}\right),$

with $\beta_{k}\in\mathbb{R}^{\underline{q}_{k}}\text{,
}\eta_{k}:\Theta\rightarrow\mathbb{R}^{\underline{d}_{k}}\text{
}g_{k}\in\mathbb{R}^{\underline{d}_{k}}$, and network change statistics
$\Delta_{i,j,k}=((g_{k}(x_{i,j}^{+})-g_{k}(x_{i,j}^{-}))^{\intercal},|x_{i,j}^{+}|-|x_{i,j}^{-}|)^{\intercal}$,
where $x_{i,j}^{+}$ and $x_{i,j}^{-}$ is network $x$ with $x_{i,j}=x_{j,i}=1$
and $x_{i,j}=x_{j,i}=0$ respectively. Likewise, MPLEs can be obtained from a
directed network, or from a weighted (undirected or directed) network based on
a binary representation of a polytomous network (Caimo & Gollini, 2022, §4.3).
Such a logistic regression likelihood a special form of composite likelihood
(Lindsay 1988; Varin et al. 2011) which assumes (dyadic) independence of the
$x_{i,j}$ observations, $\Pr_{\theta}(x_{i,j}=1\mid
X_{i,j}^{c}=x_{i,j}^{c})=\Pr_{\theta}(x_{i,j}=1)$, where $x_{i,j}^{c}$ is the
network $x$ excluding $x_{i,j}$. Both the MPLE and MLE of the ERGM are
consistent for a growing number of networks observed from the same set of
fixed nodes (Arnold & Strauss, 1991). The MPLE can be rapidly computed from
large networks (Schmid & Desmarais, 2017), if necessary, using divide-and-
conquer (Gao et al. 2022; Rosenblatt & Nadler, 2016; Minsker, 2019) or
streaming methods (Luo & Song, 2019).

###  3.3 Simulation Study: Exponential Random Graph Model (ERGM)

Here we consider a simulation study to evaluate proposed methods for the ERGM
using ABC Algorithm 1. The study involves 10 replications of $n=50$ node
undirected networks, and of $n=300$ node undirected networks simulated from
this model.

Table 2: Results of simulation study for the ERGM.
---
Mean (standard deviation) of posterior estimate over 10 replicas.
| kstar(2) | triangles | kstar(2) | triangles | degree1.5
n | 50 | 50 | 300 | 300 | 300
n${}_{\text{sim}}$ | 25 | 25 | 100 | 100 | 100
true $\theta$ | -0.20 | 0.50 | -0.20 | 0.50 | 0.80
ABC mean | -0.22 (0.03) | 0.95 (0.35) | -0.03 (0.02) | 0.49 (0.10) | -0.15 (0.12)
ABC med. | -0.22 (0.03) | 0.86 (0.33) | 0.01 (0.03) | 0.28 (0.07) | -0.85 (0.28)
ABC mode | -0.18 (0.06) | 0.57 (0.51) | 0.06 (0.07) | 1.18 (0.64) | -1.14 (0.07)
ABC MLE | -0.22 (0.06) | 0.79 (0.48) | 0.01 (0.07) | 1.61 (0.58) | -1.09 (0.09)
MCMLE | -0.20 (0.02) | 0.46 (0.16) | -0.03 (0.01) | 0.09 (0.01) | 0.08 (0.03)
MPLE | 0.15 (0.00) | 1.18 (0.08) | 0.03 (0.00) | 0.26 (0.00) | 0.18 (0.00)
ABC s.d. | 0.07 (0.01) | 0.66 (0.12) | 0.27 (0.02) | 0.79 (0.10) | 1.51 (0.21)
SE(MCMLE) | 0.02 (0.00) | 0.19 (0.01) | 0.00 (0.00) | 0.01 (0.00) | 0.02 (0.00)
Copula | d.f. | Correlation | d.f. | Correlation Matrix
d.f. & corr. | 14.55 (13.33) | -0.48 (0.12) | 20.48 (24.78) | triangles | degree1.5
|  |  | kstar(2) | -0.07 (0.08) | -0.97 (0.01)
|  |  | triangles |  | -0.09 (0.09)
MAE (MSE) of posterior estimate over 10 replicas.
| kstar(2) | triangles | kstar(2) | triangles | degree1.5
n | 50 | 50 | 300 | 300 | 300
n${}_{\text{sim}}$ | 25 | 25 | 100 | 100 | 100
ABC mean | 0.03 (0.00) | 0.47 (0.31) | 0.17 (0.03) | 0.07 (0.01) | 0.95 (0.92)
ABC med. | 0.03 (0.00) | 0.42 (0.23) | 0.21 (0.04) | 0.22 (0.05) | 1.65 (2.80)
ABC mode | 0.05 (0.00) | 0.41 (0.24) | 0.26 (0.07) | 0.72 (0.83) | 1.94 (3.77)
ABC MLE | 0.05 (0.00) | 0.43 (0.29) | 0.21 (0.05) | 1.11 (1.53) | 1.89 (3.60)
MCMLE | 0.02 (0.00) | 0.12 (0.02) | 0.17 (0.03) | 0.41 (0.16) | 0.72 (0.52)
MPLE | 0.35 (0.12) | 0.68 (0.47) | 0.23 (0.05) | 0.24 (0.06) | 0.62 (0.38)
|  |  |  |  |

Under the $n=50$ simulation condition, each of the 10 network datasets was
simulated from the ERGM (3.1) using true data-generating parameters
$\theta^{\intercal}=(-0.20,0.50)$ and defined by network sufficient statistics
being the number of 2-stars and triangles. For each simulated network dataset
$x$ analyzed by ABC Algorithm 1 using the ERGM, this model was assigned a $g$
prior
$\theta\sim\mathcal{N}_{2}(0,g(-\mathcal{H}(\widehat{\theta}_{\text{MPLE}}))^{-1})$
with $g=100$ where $\mathcal{H}(\widehat{\theta}_{\text{MPLE}})$ is the
Hessian matrix of the MPLE of $\theta$, and the summary statistics were
specified as
$s(x)=(\widehat{\beta}_{\text{MPLE},1}^{(x)},\widehat{\beta}_{\text{MPLE},2}^{(x)})^{\intercal}$,
where $\widehat{\beta}_{\text{MPLE},k}$ for $k=1,2$ are the MPLEs for these
two network sufficient statistics respectively based on the edge count $|x|$
offset with fixed coefficient $\log(1/n)$. Likewise for
$s(y)=(\widehat{\beta}_{\text{MPLE},1}^{(y)},\widehat{\beta}_{\text{MPLE},2}^{(y)})^{\intercal}$
based on the edge count $|y|$ offset with fixed coefficient
$\log(1/n_{\text{sim}})$, and on the size $n_{\text{sim}}=25$ of the network
dataset $y$ simulated in each iteration of ABC Algorithm 1.

Similarly, under the $n=300$ simulation condition, each of the 10 network
datasets was simulated from the ERGM (3.1) using true data-generating
parameters $\theta^{\intercal}=(-0.20,0.50,0.80)$ and defined by network
sufficient statistics being the number of 2-stars (kstar(2)), the number of
triangles (triangles), and the geometric weighted degree distribution
(degree1.5) with decay parameter fixed at $\alpha\equiv\log(1.5)$ (Snijders et
al. 2006, from equations 11-12, pp. 112, 126). For each simulated network
dataset $x$ analyzed by ABC Algorithm 1 using the ERGM, this model was
assigned $g$ prior
$\theta\sim\mathcal{N}_{3}(0,g(-\mathcal{H}(\widehat{\theta}_{\text{MPLE}}))^{-1})$
with $g=100,000,$ and the summary statistics were specified as
$s(x)=(\widehat{\beta}_{\text{MPLE},1}^{(x)},\widehat{\beta}_{\text{MPLE},2}^{(x)},\widehat{\beta}_{\text{MPLE},3}^{(x)})^{\intercal}$,
where $\widehat{\beta}_{\text{MPLE},k}$ for $k=1,2,3,$ are the MPLEs for these
three network sufficient statistics respectively based on the edge count $|x|$
offset with fixed coefficient $\log(1/n)$. Likewise for
$s(y)=(\widehat{\beta}_{\text{MPLE},1}^{(y)},\widehat{\beta}_{\text{MPLE},2}^{(y)},\widehat{\beta}_{\text{MPLE},3}^{(y)})^{\intercal}$
based on the edge count $|y|$ offset with fixed coefficient
$\log(1/n_{\text{sim}})$, and on the size $n_{\text{sim}}=100$ of the network
dataset $y$ simulated in each iteration of ABC Algorithm 1.

The ERGM simulation study results in Table 2 show that the MAE and MSE were
small for most ABC posterior and MLE estimators computed from networks
simulated from the model of sizes $n_{\text{sim}}<n$, which were remarkably
competitive with the MCMLE and MPLE estimators of the ERGM parameters $\theta$
each computed on the full network datasets on $n$ nodes. The ERGM results for
$n=300$ can be improved by using a smaller $g$ for the $g$-prior. The MCMLE
estimator (introduced by Snijders, 2002) is a commonly used MCMC approximation
of the MLE for the ERGM, and thus provides a standard for comparison. MCMLE is
too computationally costly and inapplicable on sufficiently large or massive
networks, because it can take a prohibitively long time to simulate a very
large network using MCMC conditionally on ERGM parameters.

###  3.4 Simulation Study: Preferential Attachment (PA) Models

We study the Price model and the NLPA model, each defined by preferential
attachment rules. Such a model, in its basic form, grows a directed or
undirected network by introducing a new node at each stage and linking the new
node to any given existing node $i$ with probability $\Pr($New node attaches
to old node $i$) $\propto k_{0}+k_{i}^{\alpha}$, with some constant parameter
$k_{0}\in\mathbb{R}$, and $k_{i}$ is the degree of node $i$ (indegree for a
directed network); and power parameter $\alpha\in(0,\infty)$.

The Price model grows a directed binary citation network one article (node) at
a time (each edge a paper citing another paper), with $\Pr($new article cites
existing article $i)\propto k_{0}+k_{i}$, with $\alpha\equiv 1$ and constant
parameter $k_{0}$ centered around 1, and $k_{i}$ is the in-degree $k_{i}$ of
existing article $i$, such that each new article cites $m$ existing articles
on average, and the number of articles a new article cites follows a
Binomial$(n_{0},p)$ distribution with success probability parameter $p$, with
$n_{0}$ the maximum outdegree of $x$ (Raynal et al. 2022). From an observed
network dataset $x$, the unknown Price model parameters to be estimated are
$(k_{0},p)$. For a undirected (binary) network $x$, the Barabasi & Albert
(1999) (BA) model has the same form except with the usual assumptions
$k_{0}=0$ and $\alpha\equiv 1$, and that $k_{i}$ is the degree of existing
node $i$.

Table 3: Results of simulation study for Price and NLPA models.
---
Mean (standard deviation) of posterior estimate over 10 replicas.
Price model | $k_{0}$ | $p$ | $k_{0}$ | $p$
n | 50 | 50 | 300 | 300
n${}_{\text{sim}}$ | 25 | 25 | 100 | 100
true $\theta$ | 1.00 | 0.02 | 1.00 | 0.02
ABC mean | 1.00 (0.00) | 0.05 (0.01) | 0.99 (0.01) | 0.05 (0.00)
ABC med. | 1.00 (0.01) | 0.05 (0.01) | 0.99 (0.01) | 0.05 (0.00)
ABC mode | 0.97 (0.05) | 0.05 (0.02) | 1.02 (0.06) | 0.05 (0.00)
MLE | 0.97 (0.05) | 0.05 (0.02) | 1.02 (0.06) | 0.05 (0.00)
ABC s.d. | 0.06 (0.00) | 0.02 (0.00) | 0.06 (0.00) | 0.01 (0.00)
Copula | d.f. | Correlation | d.f. | Correlation
d.f. & corr. | 633.32 (474.75) | 0.06 (0.05) | 568.93 (456.18) | 0.04 (0.07)
MAE (MSE) of posterior estimate over 10 replicas.
Price model | $k_{0}$ | $p$ | $k_{0}$ | $p$
n | 50 | 50 | 300 | 300
n${}_{\text{sim}}$ | 25 | 25 | 100 | 100
ABC mean | 0.00 (0.00) | 0.03 (0.00) | 0.01 (0.00) | 0.03 (0.00)
ABC med. | 0.00 (0.00) | 0.03 (0.00) | 0.01 (0.00) | 0.03 (0.00)
ABC mode | 0.05 (0.00) | 0.03 (0.00) | 0.05 (0.00) | 0.03 (0.00)
ABC MLE | 0.05 (0.00) | 0.03 (0.00) | 0.05 (0.00) | 0.03 (0.00)
Mean (standard deviation) of posterior estimate over 10 replicas.
NLPA model | $\alpha$ | $p$ | $\alpha$ | $p$
n | 50 | 50 | 300 | 300
n${}_{\text{sim}}$ | 25 | 25 | 100 | 100
true $\theta$ | 1.20 | 0.02 | 1.20 | 0.02
ABC mean | 1.24 (0.30) | 0.01 (0.00) | 1.12 (0.09) | 0.02 (0.00)
ABC med. | 1.23 (0.30) | 0.01 (0.00) | 1.12 (0.08) | 0.02 (0.00)
ABC mode | 1.35 (0.50) | 0.01 (0.00) | 1.11 (0.22) | 0.02 (0.00)
MLE | 1.35 (0.50) | 0.01 (0.00) | 1.11 (0.22) | 0.02 (0.00)
ABC s.d. | 0.51 (0.02) | 0.01 (0.00) | 0.27 (0.04) | 0.00 (0.00)
Copula | d.f. | Correlation | d.f. | Correlation
d.f. & corr. | 13.52 (9.78) | 0.21 (0.25) | 7.39 (9.64) | 0.10 (0.10)
MAE (MSE) of posterior estimate over 10 replicas.
NLPA model | $\alpha$ | $p$ | $\alpha$ | $p$
n | 50 | 50 | 300 | 300
n${}_{\text{sim}}$ | 25 | 25 | 100 | 100
ABC mean | 0.23 (0.08) | 0.01 (0.00) | 0.09 (0.01) | 0.00 (0.00)
ABC med. | 0.23 (0.08) | 0.02 (0.00) | 0.08 (0.01) | 0.00 (0.00)
ABC mode | 0.41 (0.25) | 0.01 (0.00) | 0.16 (0.05) | 0.00 (0.00)
ABC MLE | 0.41 (0.25) | 0.01 (0.00) | 0.16 (0.05) | 0.00 (0.00)
|  |  |  |

The nonlinear preferential attachment (NLPA) model (Krapivsky et al. 2000)
generalizes the BA model, using $\Pr($new node attaches to existing node
$i)=\dfrac{k_{i}^{\alpha}}{{\textstyle\sum\nolimits_{\text{Existing nodes
}j}}k_{j}^{\alpha}},$ $\alpha>0$, with $k_{i}$ the degree of existing node
$i$. Since possibly $\alpha\neq 1$, it makes sense to have the number of nodes
a new node attaches to follow a truncated Binomial$(n_{0},p)$ distribution
which supports only positive counts. For the NLPA model, the unknown
parameters to be estimated from the given network dataset are $(\alpha,p)$.
Thus from an observed network dataset $x$, the unknown model parameters to be
estimated are $(\alpha,p)$, and $n_{0}$ is chosen as the maximum degree of
$x$.

Table 3 presents results of the simulation study of ABC Algorithm 1 for the
Price and NLPA models. This study is based on 10 replications of $n=50$ node
directed networks, and 10 replications of $n=300$ node directed networks
simulated from the Price model with true data-generating parameters
$\theta^{\intercal}=(k_{0},p)=(1,0.02)$. Also, this study is based on 10
replications of $n=50$ node undirected networks, and 10 replications of
$n=300$ node undirected networks simulated from the NLPA model with true data-
generating parameters $\theta^{\intercal}=(\alpha,p)=(1.2,0.02)$. For each
network dataset $x$ simulated from the Price model and analyzed by ABC
Algorithm 1 using this model, the Price model was assigned a uniform prior
$(k_{0},p)\sim$ $\mathcal{U}_{(0.9,1.1)}\mathcal{U}_{(0,b)}$ (with $b=0.10$
under $n=50$ and $b=0.20$ for $n=300$) and using vector of summaries
$s(\cdot)$ specified by network size invariant offset-MPLEs of geometrically
weighted degree, its decay estimate, and triangle count. For each network
dataset $x$ simulated from the NLPA model and analyzed by ABC Algorithm 1
using this model, the NLPA model was assigned a uniform prior
$(\alpha,p)\sim\mathcal{U}_{(0,3)}\mathcal{U}_{(0,0.20)}$ and using vector of
summaries $s(\cdot)$ specified by the network size invariant offset-MPLEs of
density, average clustering coefficient, and diameter. The MPLE summaries
$s(x)$ of each network $x$ used offsets based on $n$ nodes ($n=50$ or $300$),
while the MPLE summaries $s(y)$ of each network $y$ of smaller size
$n_{\text{sim}}<n$ (being $n_{\text{sim}}=25$ for $n=50$, and
$n_{\text{sim}}=100$ for $n=300$) was simulated in each iteration of the ABC
algorithm used offsets based on $n_{\text{sim}}$ nodes, and based on a
Binomial$(n_{\text{sim}}-1,p)$ distribution (or the truncated binomial for
NLPA). The simulation study results of Table 3 shows that for both models the
MAE and MSE were small for the ABC posterior and MLE estimators computed from
networks simulated from the model of sizes $n_{\text{sim}}<n$.

###  3.5 Simulation Study: Duplication-Divergence Models

We now study the DMC model and the DMR model. Each model grows an undirected
network of $n$ nodes by starting with a seed network of $n^{\prime}\geq 1$
nodes, and then repeats the following steps until requisite number of nodes
($n=n^{\prime}$) is reached: Add new node ($n^{\prime}\equiv n^{\prime}+1$),
then add edge between the new node and each neighbor of a randomly chosen
existing node, $k\sim\mathcal{U}_{\\{1,n^{\prime}-1\\}}$; then under DMC: for
each neighbor of the chosen node, randomly select either the edge between the
chosen node and the neighbor, or the edge between the new node and the
neighbor, and remove that edge with probability $q_{\operatorname{mod}}$; and
then add an edge between the chosen node and the new node with probability
$q_{\text{con}}$; and instead, under DMR: each edge connected to the new node
is removed independently with probability $q_{\text{del}}$, and an edge
between any existing node and the new node is added with probability
$q_{\text{new}}/(n^{\prime}-1)$.

Table 4: Results of simulation study for the DMC and DMR models.
---
Mean (standard deviation) of posterior estimate over 10 replicas.
DMC model | $q_{\operatorname{mod}}$ | $q_{\text{con}}$ | $q_{\operatorname{mod}}$ | $q_{\text{con}}$
n | 50 | 50 | 300 | 300
n${}_{\text{sim}}$ | 25 | 25 | 100 | 100
true $\theta$ | 0.20 | 0.10 | 0.20 | 0.10
ABC mean | 0.24 (0.02) | 0.16 (0.05) | 0.22 (0.03) | 0.16 (0.05)
ABC med. | 0.24 (0.03) | 0.15 (0.05) | 0.22 (0.04) | 0.14 (0.04)
ABC mode | 0.23 (0.06) | 0.16 (0.09) | 0.26 (0.04) | 0.13 (0.09)
MLE | 0.23 (0.06) | 0.16 (0.09) | 0.26 (0.04) | 0.13 (0.09)
ABC s.d. | 0.05 (0.01) | 0.07 (0.02) | 0.04 (0.01) | 0.07 (0.04)
Copula | d.f. | Correlation | d.f. | Correlation
d.f. & corr. | 360.66 (448.23) | 0.09 (0.24) | 131.75 (289.23) | 0.21 (0.30)
MAE (MSE) of posterior estimate over 10 replicas.
DMC model | $q_{\operatorname{mod}}$ | $q_{\text{con}}$ | $q_{\operatorname{mod}}$ | $q_{\text{con}}$
n | 50 | 50 | 300 | 300
n${}_{\text{sim}}$ | 25 | 25 | 100 | 100
ABC mean | 0.04 (0.00) | 0.06 (0.01) | 0.03 (0.00) | 0.06 (0.01)
ABC med. | 0.04 (0.00) | 0.05 (0.01) | 0.03 (0.00) | 0.05 (0.00)
ABC mode | 0.04 (0.00) | 0.08 (0.01) | 0.06 (0.00) | 0.06 (0.01)
ABC MLE | 0.04 (0.00) | 0.08 (0.01) | 0.06 (0.00) | 0.06 (0.01)
Mean (standard deviation) of posterior estimate over 10 replicas.
DMR model | $q_{\text{del}}$ | $q_{\text{new}}$ | $q_{\text{del}}$ | $q_{\text{new}}$
n | 50 | 50 | 300 | 300
n${}_{\text{sim}}$ | 25 | 25 | 100 | 100
true $\theta$ |  0.20 | 0.10 |  0.20 | 0.10
ABC mean | 0.24 (0.03) | 0.24 (0.19) | 0.25 (0.01) | 0.13 (0.06)
ABC med. | 0.24 (0.03) | 0.23 (0.21) | 0.25 (0.02) | 0.12 (0.06)
ABC mode | 0.23 (0.06) | 0.28 (0.20) | 0.26 (0.05) | 0.12 (0.09)
MLE | 0.23 (0.06) | 0.28 (0.20) | 0.26 (0.05) | 0.12 (0.09)
ABC s.d. | 0.05 (0.01) | 0.12 (0.04) | 0.06 (0.00) | 0.07 (0.02)
Copula | d.f. | Correlation | d.f. | Correlation
d.f. & corr. | 426.42 (425.06) | -0.05 (0.19) | 354.68 (451.91) | 0.02 (0.17)
MAE (MSE) of posterior estimate over 10 replicas.
DMR model | $q_{\text{del}}$ | $q_{\text{new}}$ | $q_{\text{del}}$ | $q_{\text{new}}$
n | 50 | 50 | 300 | 300
n${}_{\text{sim}}$ | 25 | 25 | 100 | 100
ABC mean | 0.05 (0.00) | 0.14 (0.05) | 0.05 (0.00) | 0.04 (0.00)
ABC med. | 0.05 (0.00) | 0.13 (0.06) | 0.05 (0.00) | 0.05 (0.00)
ABC mode | 0.05 (0.00) | 0.19 (0.07) | 0.06 (0.01) | 0.07 (0.01)
ABC MLE | 0.05 (0.00) | 0.19 (0.07) | 0.06 (0.01) | 0.07 (0.01)
|  |  |  |

Table 4 presents results of the simulation study of ABC Algorithm 1 for the
DMC and DMR models. This study is based on 10 replications of $n=50$ node
undirected networks, and 10 replications of $n=300$ node undirected networks
simulated from the DMC model with true data-generating parameters
$\theta^{\intercal}=(q_{\operatorname{mod}},q_{\text{con}})=(0.20,0.10)$.
Also, this study is based on 10 replications of $n=50$ node undirected
networks, and 10 replications of $n=300$ node undirected networks simulated
from the DMR model with true data-generating parameters
$\theta^{\intercal}=(q_{\text{del}},q_{\text{new}})=(0.20,0.10)$. For each
network dataset $x$ simulated from the DMC model (DMR model, respectively) and
analyzed by ABC Algorithm 1 using this model, the DMC model (DMR model, resp.)
was assigned a uniform prior
$(q_{\operatorname{mod}},q_{\text{con}})\sim\mathcal{U}_{(0.15,0.35)}\mathcal{U}_{(0,1)}$
(uniform prior
$(q_{\text{del}},q_{\text{new}})\sim\mathcal{U}_{(0.15,0.35)}\mathcal{U}_{(0,1)}$,
resp.) and using vector of summaries $s(\cdot)$ specified by network size
invariant offset-MPLEs of undirected network summary statistics of mean degree
and triangles (of local and global average clustering coefficients, and degree
assortativity, resp.). The MPLE summaries $s(x)$ of each network $x$ used
offsets based on $n$ nodes ($n=50$ or $300$), while the MPLE summaries $s(y)$
of each network $y$ of smaller size $n_{\text{sim}}<n$ (being
$n_{\text{sim}}=25$ for $n=50$, and $n_{\text{sim}}=100$ for $n=300$) was
simulated in each iteration of the ABC algorithm used offsets based on
$n_{\text{sim}}$ nodes. Again, the simulation study results of Table 4 shows
that for both models the MAE and MSE were small for the ABC posterior and MLE
estimators computed from networks simulated from the model of sizes
$n_{\text{sim}}<n$.

###  3.6 Real Citation Network Analysis

With ABC Algorithm 1 validated by simulation studies for several network
models, we turn to applications of this algorithm to analyze real-life large
network datasets. We now analyze a large binary directed citation HepPh
network dataset, a Arxiv High Energy Physics paper citation network of
$n=28,093$ papers (nodes) and $3,148,447$ citations (directed edges) among
them (Rossi & Ahmed, 2015). (Dataset obtained from
https://networkrepository.com/cit.php). This dataset is analyzed by the ERGM
model and the Price model, each model using a vector of summary statistics
$s(\cdot)$ being the network size invariant offset-MPLEs (resp.) of
geometrically weighted in-degree distribution (gwidegree), its decay (degree
weighting) parameter estimate (gwidegree.decay), and triangle count. The ERGM,
based on network sufficient statistics gwidegree, gwidegree.decay, and
triangle, was assigned a trivariate normal prior $\theta\sim$
$\mathcal{N}(0,10I_{3})$. The Price model was assigned priors $k_{0}\sim$
$\mathcal{U}_{(0.9,1.1)}$ and $p\sim\mathcal{U}_{(0,0.20)}$.

It turned out that the MPLEs for the three network summaries (resp.) could not
be computed directly on this large network, as the computations were
prohibitively long. Therefore, these 3 MPLEs were estimated by their geometric
median over 100,000 subsamples of the 28,093 papers (nodes), each subsample
inducing a subgraph of size $\lfloor\sqrt{28,093}\rfloor=167$, providing an
outlier-robust, divide-and-conquer (DAC) subsampling MPLE estimator (Minsker
2019). This resulted in the following MPLEs (resp.): for gwidegree,
-0.8418461; for gwidegree.decay, 0.2365466; for triangle: 1.6904300. Recall
from §3.2 that the MPLE is the MLE for logistic regression assuming dyadic
independence of the edge observations, which are both consistent over a
growing number of networks observed from the same set of fixed nodes. These
three DAC MPLEs then specified the vector of summary statistics $s(x)$ for the
real observed citation HepPh network dataset, analyzed by each of the ERGM and
Price model using ABC Algorithm 1, with each model simulating a network of
size $n_{\text{sim}}=$ $167$ in each algorithm iteration.

Table 5 presents the posterior estimates and MLEs obtained from this ABC
algorithm. Running Step 2 of this algorithm, for model selection, output a
posterior probability estimate of $0.99$ for the Price model, compared to the
ERGM.

Table 5: ABC posterior estimates for the Price model and the ERGM.
---
| Price model |  | ERGM |
| $k_{0}$ | $p$ |  $\theta$ gwidegree |  $\theta$ decay |  $\theta$ triangles
ABC mean | 1.01 | 0.01 | -4.46 | 3.34 | -1.07
ABC median | 1.02 | 0.01 | -4.96 | 3.66 | -0.58
ABC mode | 1.08 | 0.01 | -1.76 | -0.80 | -4.59
MLE | 1.08 | 0.01 | -4.82 | 6.00 | -0.58
ABC s.d. | 0.05 | 0.00 | 2.41 | 2.62 | 2.06
ABC IQR |  (0.97, 1.06) |  (0.01, 0.01) | (-6.09, -2.81) | (2.23, 4.96) | (-1.17, -0.23)
ABC 95% |  (0.91, 1.09) |  (0.00, 0.02) | (-9.12, 0.73) | (-2.08, 6.80) | (-6.09, 3.22)
Copula | d.f. | Correlation | d.f. | Correlation matrix
d.f. & corr. | 50.42 | -0.20 | 6.31 |  $\theta$ decay |  $\theta$ triangles
|  |  |  $\theta$ gwidegree | -0.60 | 0.30
|  |  |  $\theta$ decay |  | -0.28

###  3.7 Real Multilayer Network Analysis

The multilayer, weighted, and directed, BostonBomb2013 Twitter network dataset
consists of $n\ =4,377,184$ persons (nodes), and $9,480,331$ count-weighted
edges, across $L=3$ layers: retweets, mentions and replies, occurring between
dates 4-15-2013 and 4-22-2013 of the Boston Bombing Attacks of 2013 (De
Domenico & Altmann 2020). (Dataset obtained from
https://manliodedomenico.com/data.php). Using the ABC Algorithm 1, the
BostonBomb2013 network dataset was analyzed by a three-Layer ERGM with Poisson
reference measure $b$,

$f(X=x\mid\theta)={\displaystyle\prod\limits_{l=1}^{L=3}}f(X_{l}=x_{l}\mid\theta)=\dfrac{b(x_{l})\exp\left\\{{\textstyle\sum\nolimits_{l=1}^{L=3}}\theta_{l}^{\intercal}h_{l}(x_{l})\right\\}}{{\textstyle\prod\limits_{l=1}^{L=3}}{\textstyle\sum\limits_{x^{\prime}\in\mathcal{X}}}b(x^{\prime})\exp\left\\{{\textstyle\sum\nolimits_{l=1}^{L=3}}\theta_{l}^{\intercal}h_{l}(x^{\prime})\right\\}},$

which was assigned a multivariate normal prior prior $\theta\sim$
$\mathcal{N}(0,10I_{18})$ for 18 parameters $\theta$ of six network sufficient
statistics (vector $h_{l}$) specified for each of the 3 layers, namely: (1)
Indicator of edge weight equal to 1 in the given layer (equalto(1)); (2)
Indicator of edge weight greater than 1 in the given layer; (3) minimum
mutuality (mutual.min) which is analogous to mutuality for binary networks;
(4) triad closure represented by transitive weights with twopath = min,
combine = sum, and affect = min (transitiveweights.min.sum.min), which is
analogous to triangle count from a binary network; (5) triad closure
represented by transitive weights with twopath = min, combine = max, and
affect = min (transitiveweights.min.max.min), which is analogous to a count of
transitive ties from a binary network; and the term CMP, which specifies the
reference measure $b$ by the Conway-Maxwell-Poisson (CMP) distribution for the
count-weighted edges, with corresponding coefficient parameter
$\theta_{\text{CMP}}$, which controls the degree of dispersion relative to
Poisson distribution. In particular, $\theta_{\text{CMP}}=$ 0 defines the
Poisson distribution; $\ \theta_{\text{CMP}}$ $<0$ more underdispersed
distribution; $\theta_{\text{CMP}}$ $\rightarrow$ $-\infty$ approaches
Bernoulli distribution; $\theta_{\text{CMP}}$ $>0$ more overdispersed
distribution; $\theta_{\text{CMP}}=1$ geometric distribution (most
overdispersed). More details about these network statistics of valued networks
and the CMP distribution for ERGM are provided elsewhere(Krivitsky, 2012).

Table 6: Summary statistics $s(x)$ of the BostonBomb2013 network dataset.
---
Summary Statistic | Layer 1 | Layer 2 | Layer 3
meanIndegree | 1.35253692785 | 0.6141674647 | 0.1590769317
varIndegree | 8474.74052782986 | 1804.0675198631 | 31.5665679809
meanOutdegree | 1.35253692785 | 0.6141674647 | 0.1590769317
varOutdegree | 7.61861428910 | 4.5875615155 | 0.8667280838
wClusteringCoef | 0.00008622621 | 0.0004401573 | 0.0001867074
assortativityDegree | -0.06229008973 | -0.0233535797 | -0.0165706996
reciprocity | 0.00379811641 | 0.0379874337 | 0.0625188215

Table 7: ABC posterior estimates of the 3-layer valued ERGM parameters
---
from the BostonBomb2013 network dataset (mode and MLE identical).
Layer 1 | Mean | Med. | MLE | S.D. | 2.5% | 25% | 75% | 97.5%
equalto.1.pm.0 | -5.65 | -5.42 | -2.94 | 1.38 | -9.07 | -6.30 | -4.69 | -3.59
greaterthan.1 | -2.85 | -3.18 | 6.32 | 2.74 | -7.72 | -4.77 | -1.06 | 3.06
mutual.min | -3.85 | -3.63 | 1.20 | 2.20 | -8.72 | -5.27 | -2.42 | 0.14
tw.min.sum.min | 0.05 | 0.18 | -0.22 | 3.25 | -6.71 | -2.30 | 2.38 | 5.74
tw.min.max.min | -1.88 | -1.90 | 3.79 | 2.89 | -7.29 | -3.84 | 0.06 | 3.78
CMP | -2.15 | -2.24 | 3.77 | 2.66 | -7.33 | -3.89 | -0.38 | 3.04
Layer 2 | Mean | 50% | MLE | S.D. | 2.5% | 25% | 75% | 97.5%
equalto.1.pm.0 | -6.39 | -6.33 | -3.12 | 1.44 | -9.76 | -7.26 | -5.39 | -3.91
greaterthan.1 | -2.76 | -2.66 | 8.51 | 2.80 | -7.83 | -4.71 | -1.02 | 3.13
mutual.min | 1.70 | 1.59 | -3.52 | 1.25 | -0.32 | 0.78 | 2.57 | 3.83
tw.min.sum.min | 0.23 | 0.27 | 1.38 | 2.98 | -5.81 | -1.64 | 2.33 | 5.87
tw.min.max.min | -2.48 | -2.45 | 5.79 | 2.79 | -8.02 | -4.36 | -0.64 | 3.30
CMP | -2.04 | -2.10 | 3.24 | 2.93 | -7.55 | -3.95 | -0.01 | 3.84
Layer 3 | Mean | 50% | MLE | S.D. | 2.5% | 25% | 75% | 97.5%
equalto.1.pm.0 | -5.67 | -5.62 | -3.04 | 1.25 | -8.44 | -6.52 | -4.78 | -3.71
greaterthan.1 | -3.22 | -3.17 | 6.77 | 2.76 | -8.22 | -5.20 | -1.45 | 2.37
mutual.min | 2.62 | 2.41 | -2.55 | 1.21 | 0.65 | 1.86 | 3.33 | 5.16
tw.min.sum.min | -1.14 | -1.23 | 2.96 | 2.82 | -6.86 | -2.79 | 0.68 | 4.23
tw.min.max.min | -2.15 | -2.28 | 8.49 | 2.85 | -7.68 | -4.03 | -0.33 | 3.92
CMP | -1.34 | -1.20 | 5.61 | 2.72 | -6.82 | -3.20 | 0.47 | 4.11
Notes: tw.min.sum.min is transitiveweights.min.sum.min.
tw.min.max.min refers to transitiveweights.min.max.min.

These 18 ERGM sufficient statistics could not be easily computed on the entire
BostonBomb2013 network dataset, due to its massive size. Therefore, to
approximate these sufficient statistics, the vector of summary statistics
$s(x)$ was computed by the following 21 (mostly) UMP summary statistics, each
of which were computable on the full network, for each of the 3 network
layers. The values of the summaries $s(x)$ of the BostonBomb2013 network
dataset are shown in Table 6. The ABC Algorithm 1 was run on this model based
on these network data summaries $s(x)$ and based on summaries $s(y)$ of each
network $y$ of $n_{\text{sim}}=500$ nodes simulated from the ERGM
conditionally on given proposed model parameters, in each iteration of this
algorithm. The posterior distribution and MLE estimates delivered by the
algorithm are shown in Table 7.

###  3.8 Real Social Network Analysis

Finally, we apply ABC Algorithm 1 to analyze the massive friendster social
network of $n=65,608,366$ persons (nodes) and 1,806,067,135 undirected edges
(Yang & Leskovec, 2012) (data from https://snap.stanford.edu/data/com-
Friendster.html). The NLPA model, assigned uniform priors
$(\alpha,p)\sim\mathcal{U}_{(0,3)}\mathcal{U}_{(0,0.20)}$, was used to analyze
this network dataset $x$, with summary statistics $s(x)=$ (Average clustering
coefficient, Diameter, Density) $=(0.1623,32,8.39161\times 10^{-7})$, and with
these network summaries also computed from a network of size
$n_{\text{sim}}=1,000$ simulated from the model (given a set of proposed
parameters) in each of the $N=10,000$ iteration of the ABC algorithm. The
algorithm routinely delivered posterior distribution and MLE estimates shown
in Table 8, despite the sheer size of the dataset.

Table 8: Posterior estimates of the NLPA model
---
for the friendster network dataset.
| $\alpha$ | $p$
ABC mean | 1.21 | 0.001
ABC median | 1.25 | 0.001
ABC mode | 0.24 | 0.0002
MLE | 0.24 | 0.0002
ABC s.d. | 0.28 | 0.002
ABC IQR |  (1.18, 1.31) |  (0.0003, 0.0012)
ABC 95% |  (0.37, 2.05) |  (0.0001, 0.0046)
Copula | d.f. 28.91 | Corr. 0.32

## Acknowledgements

The author thanks the journal’s Editorial Board for considering the review of
this manuscript. The code used to generate the results of this manuscript can
be made available upon request of the author.

## References

Abraham, C., Biau, G., and Cadre, B. (2003). Simple estimation of the mode of
a multivariate density. Canadian Journal of Statistics, 31: 23-34.

Abraham, C., Biau, G., and Cadre, B. (2004). ”On the asymptotic properties of
a simple estimate of the mode.” ESAIM: Probability and Statistics, 8: 1-11.

Acar, E., Craiu, R., and Yao, F. (2013). ”Statistical testing of covariate
effects in conditional copula models.” Electronic Journal of Statistics, 7:
2822-2850.

Arnold, B. and Strauss, D. (1991). ”Pseudolikelihood estimation: Some
examples.” Sankhya: The Indian Journal of Statistics, Series B, 53: 233-243.

Baillien, J., Gijbels, I., and Verhasselt, A. (2022). ”Estimation in copula
models with two-piece skewed margins using the inference for margins method.”
Econometrics and Statistics.

Barabasi, A. and Albert, R. (1999). ”Emergence of scaling in random networks.”
Science, 286: 509-512.

Barrientos, A. and Peña, V. (2020). ”Bayesian bootstraps for massive data.”
Bayesian Analysis, 15: 363-388.

Besag, J. (1974). ”Spatial interaction and the statistical analysis of lattice
systems (with discussion).” Journal of the Royal Statistical Society, Series
B, 36: 192-236.

Bickel, P., Klaassen, C., Ritov, Y., and Wellner, J. (1993). Efficient and
Adaptive Estimation for Semiparametric Models. Baltimore, MD: Johns Hopkins
University Press.

Breiman, L., Friedman, J., Stone, C., and Olshen, R. (1984). Classification
and Regression Trees. Boca Raton, FL: CRC Press.

Bródka, P. and Kazienko, P. (2018). ”Multilayer Social Networks.” In Alhajj,
R. and Rokne, J. (eds.), Encyclopedia of Social Network Analysis and Mining,
1408-1422. New York: Springer.

Caimo, A. and Gollini, I. (2022). ”Recent advances in exponential random graph
modelling.” Mathematical Proceedings of the Royal Irish Academy, 123: 1-12.

Ćevid, D., Michel, L., Näf, J., Bühlmann, P., and Meinshausen, N. (2022).
”Distributional random forests: Heterogeneity adjustment and multivariate
distributional regression.” Journal of Machine Learning Research, 23:
14987-15065.

Chacón, J. (2020). ”The modal age of statistics.” International Statistical
Review, 88: 122-141.

Chen, Y. and Gutmann, M. (2019). ”Adaptive Gaussian Copula ABC.” In Chaudhuri,
K. and Sugiyama, M. (eds.), Proceedings of the Twenty-Second International
Conference on Artificial Intelligence and Statistics, volume 89 of Proceedings
of Machine Learning Research, 1584-1592. Naha, Okinawa, Japan: PMLR.

Craiu, R. and Levi, E. (2022). ”Approximate methods for Bayesian computation.”
Annual Review of Statistics and Its Application, 10: 379-399.

Cranmer, K., Brehmer, J., and Louppe, G. (2020). ”The frontier of simulation-
based inference.” Proceedings of the National Academy of Sciences, 117:
30055-30062.

Csárdi, G., Nepusz, T., Traag, V., Horvt, S., Zanini, F., Noom, D., and
Müller, K. (2024). igraph: Network Analysis and Visualization in R. R package
version 1.5.1. URL https://CRAN.R-project.org/package=igraph

Dasgupta, S. and Kpotufe, S. (2014). ”Optimal rates for k-NN density and mode
estimation.” In Ghahramani, Z., Welling, M., Cortes, C., Lawrence, N., and
Weinberger, K. (eds.), Advances in Neural Information Processing Systems,
volume 27, 1-9. Curran Associates, Inc.

Daul, S., De Giorgi, E., Lindskog, F., and McNeil, A. (2003). ”The grouped
t-copula with an application to credit risk.” SSRN, 1-7.

De Domenico, M. and Altmann, E. (2020). ”Unraveling the origin of social
bursts in collective attention.” Scientific Reports, 10: 4629.

Dean, T., Singh, S., Jasra, A., and Peters, G. (2014). ”Parameter estimation
for hidden Markov models with intractable likelihoods.” Scandinavian Journal
of Statistics, 41: 970-987.

Dellaportas, P. and Tsionas, M. (2019). ”Importance sampling from posterior
distributions using copula-like approximations.” Journal of Econometrics, 210:
45-57.

Demarta, S. and McNeil, A. (2005). ”The t copula and related copulas.”
International Statistical Review, 73: 111-129.

Denuit, M. and Lambert, P. (2005). ”Constraints on concordance measures in
bivariate discrete data.” Journal of Multivariate Analysis, 93: 40-57.

Devroye, L. (1979). ”Recursive estimation of the mode of a multivariate
density.” Canadian Journal of Statistics, 7: 159-167.

Drovandi, C., Nott, D., and Frazier, D. (2023). ”Improving the accuracy of
marginal approximations in likelihood-free inference via localization.”
Journal of Computational and Graphical Statistics, 1-11.

Faisal, F., Newaz, K., Chaney, J., Li, J., Emrich, S., Clark, P., and
Milenković, T. (2017). ”GRAFENE: Graphlet-based alignment-free network
approach integrates 3D structural and sequence (residue order) data to improve
protein structural comparison.” Scientific Reports, 7: 1-15.

Fang, H., Fang, K., and Kotz, S. (2002). ”The meta-elliptical distributions
with given marginals.” Journal of Multivariate Analysis, 82: 1-16.

Frank, O. and Strauss, D. (1986). ”Markov graphs.” Journal of the American
Statistical Association, 81: 832-842.

Gao, Y., Liu, W., Wang, H., Wang, X., Yan, Y., and Zhang, R. (2022). ”A review
of distributed statistical inference.” Statistical Theory and Related Fields,
6: 89-99.

Geenens, G. (2020). ”Copula modeling for discrete random vectors.” Dependence
Modeling, 8: 417-440.

Genest, C., Ghoudi, K., and Rivest, L.-P. (1995). ”A semiparametric estimation
procedure of dependence parameters in multivariate families of distributions.”
Biometrika, 82: 543-552.

Ghafouri, S. and Khasteh, S. (2020). ”A survey on exponential random graph
models: An application perspective.” PeerJ Computer Science, 6: e269.

Goldenberg, A., Zheng, A., Fienberg, S., and Airoldi, E. (2010). ”A survey of
statistical network models.” Foundations and Trends in Machine Learning, 2:
129-233.

Gourieroux, C., Monfort, A., and Renault, E. (1993). ”Indirect inference.”
Journal of Applied Econometrics, 8: S85-S118.

Grazian, C. and Fan, Y. (2020). ”A review of approximate Bayesian computation
methods via density estimation: Inference for simulator-models.” WIREs
Computational Statistics, 12: e1486.

Gretton, A., Borgwardt, K., Rasch, M., Scholkopf, B., and Smola, A. (2007). ”A
kernel method for the two-sample-problem.” In Platt, J., Koller, D., Singer,
Y., and Roweis, S. (eds.), Advances in Neural Information Processing Systems,
volume 19, 513-520. The MIT Press.

Gutmann, M. and Corander, J. (2016). ”Bayesian optimization for likelihood-
free inference of simulator-based statistical models.” Journal of Machine
Learning Research, 17: 1-47.

Hammoud, Z. and Kramer, F. (2020). ”Multilayer networks: Aspects,
implementations, and application in biomedicine.” Big Data Analytics, 5: 1-18.

Handcock, M., Hunter, D., Butts, C., Goodreau, A., Krivitsky, P., and Morris,
M. (2023). ergm: Fit, Simulate and Diagnose Exponential-Family Models for
Networks. The Statnet Project. R package version 4.5.0.

Hanneke, S., Fu, W., and Xing, E. (2010). ”Discrete temporal models of social
networks.” Electronic Journal of Statistics, 4: 585-605.

Harris, J. (2013). An Introduction to Exponential Random Graph Modeling.
Thousand Oaks, CA: Sage.

Hintz, E., Hofert, M., and Lemieux, C. (2022a). ”Multivariate normal variance
mixtures in R: The R package nvmix.” Journal of Statistical Software, 102:
1-31.

Hintz, E., Hofert, M., and Lemieux, C. (2022b). ”Computational challenges of t
and related copulas.” Journal of Data Science, 20: 95-110.

Horvat, E. and Zweig, K. (2018). ”Multiplex Networks.” In Alhajj, R. and
Rokne, J. (eds.), Encyclopedia of Social Network Analysis and Mining,
1430-1434. New York: Springer.

Hsu, C. and Wu, T. (2013). ”Efficient estimation of the mode of continuous
multivariate data.” Computational Statistics and Data Analysis, 63: 148-159.

Hutson, A., Wilding, G., Mashtare, T., and Vexler, A. (2015). ”Measures of
biomarker dependence using a copula-based multivariate epsilon-skew-normal
family of distributions.” Journal of Applied Statistics, 42: 2734-2753.

Kajihara, T., Kanagawa, M., Yamazaki, K., and Fukumizu, K. (2018). ”Kernel
recursive ABC: Point estimation with intractable likelihood.” In Dy, J. and
Krause, A. (eds.), Proceedings of the 35th International Conference on Machine
Learning, volume 80 of Proceedings of Machine Learning Research, 2400-2409.
Stockholm, Sweden: PMLR.

Karabatsos, G. (2023). ”Approximate Bayesian computation using asymptotically
normal point estimates.” Computational Statistics, 38: 531-568.

Karabatsos, G. and Leisen, F. (2018). ”An approximate likelihood perspective
on ABC methods.” Statistics Surveys, 12: 66-104.

Kinsley, A., Rossi, G., Silk, M., and VanderWaal, K. (2020). ”Multilayer and
multiplex networks: An introduction to their use in veterinary epidemiology.”
Frontiers in Vetinary Science, 7: 596.

Kleijn, B. and van der Vaart, A. (2012). ”The Bernstein-Von-Mises theorem
under misspecification.” Electronic Journal of Statistics, 6: 354-381.

Klein, N., Stanley Smith, M., Nott, D., and Chrisholm, R. (2024). ”Regression
copulas for multivariate responses.” arXiv, e-print: 2401.11804.

Klemm, K. and Eguiluz, V. (2002). ”Highly clustered scale-free networks.”
Physical Review E, 65: 036123.

Kosmidis, I. and Karlis, D. (2016). ”Model-based clustering using copulas with
applications.” Statistics and Computing, 26: 1079-1099.

Krapivsky, P., Redner, S., and Leyvraz, F. (2000). ”Connectivity of Growing
Random Networks.” Physical Review Letters, 85: 4629-4632.

Kretzschmar, M. and Morris, M. (1996). ”Measures of concurrency in networks
and the spread of infectious disease.” Mathematical Biosciences, 133: 165-195.

Krivitsky, P. (2012). ”Exponential-family random graph models for valued
networks.” Electronic Journal of Statistics, 6: 1100-1128.

Krivitsky, P.(2022). ergm.count: Fit, simulate and diagnose exponential-family
models for networks with count edges. R package version 4.1.1.

Krivitsky, P. and Handcock, M. (2014). ”A separable model for dynamic
networks.” Journal of the Royal Statistical Society, Series B, 76: 29-46.

Krivitsky, P., Handcock, M., and Morris, M. (2011). ”Adjusting for network
size and composition effects in exponential-family random graph models.”
Statistical Methodology, 8: 319-339.

Krivitsky, P., Koehly, L., and Marcum, C. (2020). ”Exponential-family random
graph models for multi-layer networks.” Psychometrika, 85: 630-659.

Kumpula, J., Onnela, J., Saramäki, J., Kaski, K., and Kertész, J. (2007).
”Emergence of communities in weighted networks.” Physical Review Letters, 99:
228701.

Lange, K., Little, R., and Taylor, J. (1989). ”Robust statistical modeling
using the t distribution.” Journal of the American Statistical Association,
84: 881-896.

Lee, J., Li, G., and Wilson, J. (2020). ”Varying-coefficient models for
dynamic networks.” Computational Statistics and Data Analysis, 152: 107052.

Li, J., Nott, D., Fan, Y., and Sisson, S. (2017). ”Extending Approximate
Bayesian Computation methods to high dimensions via a Gaussian copula model.”
Computational Statistics and Data Analysis, 106: 77-89.

Lin, Y. and Jeon, Y. (2006). ”Random forests and adaptive nearest neighbors.”
Journal of the American Statistical Association, 101: 578-590.

Lindsay, B. (1988). ”Composite likelihood methods.” Contemporary Mathematics,
80: 221-239.

Loyal, J. and Chen, Y. (2020). ”Statistical network analysis: A review with
applications to the coronavirus disease 2019 pandemic.” International
Statistical Review, 88: 419-440.

Luo, L. and Song, P. (2019). ”Renewable estimation and incremental inference
in generalized linear models with streaming data sets.” Journal of the Royal
Statistical Society Series B, 82: 69-97.

Lusher, D., Koskinen, J., and Robins, G. (2013). Exponential Random Graph
Models for Social Networks: Theory, Methods, and Applications. U.K.: Cambridge
University Press.

Lyddon, S., Holmes, C., and Walker, S. (2019). ”General Bayesian updating and
the loss-likelihood bootstrap.” Biometrika, 106: 465-478.

Madsen, L. and Fang, Y. (2011). ”Joint regression analysis for discrete
longitudinal data.” Biometrics, 67: 1171-1176.

Marin, J., Pudlo, P., Estoup, A., and Robert, C. (2018). ”Likelihood-free
model choice.” In Sisson, S., Fan, Y., and Beaumont, M. (eds.), Handbook of
Approximate Bayesian Computation. Boca Raton, FL: CRC Press.

Martin, G., Frazier, D., and Robert, C. (2023). ”Approximating Bayes in the
21st Century.” Statistical Science, 1-26.

McFadden, D. (1989). ”A method of simulated moments for estimation of discrete
response models without numerical integration.” Econometrica, 57: 995-1026.

Michel, L. and Cevid, D. (2021). drf: Distributional Random Forests. R package
version 1.1.0.

Minsker, S. (2019). ”Distributed statistical estimation and rates of
convergence in normal approximation.” Electronic Journal of Statistics, 13:
5213-5252.

Morris, M., Goodreau, S., and Moody, J. (2007). ”Sexual networks, concurrency,
and STD/HIV.” In Holmes, K., Sparling, P., and Stamm, W. (eds.), Sexually
Transmitted Diseases, 109-126. New York: McGraw-Hill International Book Co.

Morris, M. and Kretzschmar, M. (1997). ”Concurrent partnerships and the spread
of HIV.” AIDS, 11: 641-648.

Nelsen, R. (2006). An Introduction to Copulas. New York: Springer.

Newman, M. (2003). ”The structure and function of complex networks.” SIAM
Review, 45: 167-256.

Newman, M. (2018). Networks (2nd Ed.). U.K.: Oxford University Press.

Newton, M., Polson, N., and Xu, J. (2021). ”Weighted Bayesian bootstrap for
scalable posterior distributions.” Canadian Journal of Statistics, 49:
421-437.

Nott, D., Ong, V., Fan, Y., and Sisson, S. (2018). ”High-dimensional ABC.” In
Sisson, S., Fan, Y., and Beaumont, M. (eds.), Handbook of Approximate Bayesian
Computation, 106-117. Boca Raton, FL: CRC Press.

Okhrin, O., Ristig, A., and Xu, Y. (2017). ”Copulae in High Dimensions: An
Introduction.” In Härdle, W., Chen, C., and Overbeck, L. (eds.), Applied
Quantitative Finance, 247-277. Berlin, Heidelberg: Springer Berlin Heidelberg.

Palombi, L., Bernava, G., Nucita, A., Giglio, P., Liotta, G., Nielsen-Saines,
K., Orlando, S., Mancinelli, S., Buonomo, E., Scarcella, P., Guidotti, G.,
Ceffa, S., Haswell, J., Zimba, I., Magid, N., and Marazzi, M. (2012).
”Predicting trends in HIV-1 sexual transmission in sub-Saharan Africa through
the drug resource enhancement against AIDS and malnutrition model:
Antiretrovirals for reduction of population infectivity, incidence and
prevalence at the district level.” Clinical Infectious Diseases, 55: 268-275.

Pastor-Satorras, R., Smith, E., and Solé, R. (2003). ”Evolving protein
interaction networks through gene duplication.” Journal of Theoretical
Biology, 222: 199-210.

Pastor-Satorras, R. and Vespignani, A. (2004). Evolution and Structure of the
Internet: A Statistical Physics Approach. U.K.: Cambridge University Press.

Pesonen, H., Simola, U., Kohn-Luque, A., Vuollekoski, H., Lai, X., Frigessi,
A., Kaski, S., Frazier, D., Maneesoonthorn, W., Martin, G., and Corander, J.
(2023). ”ABC of the future.” International Statistical Review, 91: 243-268.

Picchini, U. and Anderson, R. (2017). ”Approximate maximum likelihood
estimation using data-cloning ABC.” Computational Statistics and Data
Analysis, 105: 166-183.

Picchini, U., Simola, U., and Corander, J. (2022). ”Sequentially guided MCMC
proposals for synthetic likelihoods and correlated synthetic likelihoods.”
Bayesian Analysis, 18: 1099-1129.

Pitt, M., Chan, D., and Kohn, R. (2006). ”Efficient Bayesian inference for
Gaussian copula regression models.” Biometrika, 93: 537-554.

Price, D. (1965). ”Networks of scientific papers.” Science, 149: 510-515.

Price, D. (1976). ”A general theory of bibliometric and other cumulative
advantage processes.” Journal of the American Society for Information Science,
27: 292-306.

Procopio, A., Cesarelli, G., Donisi, L., Merola, A., Amato, F., and Cosentino,
C. (2023). ”Combined mechanistic modeling and machine-learning approaches in
systems biology - A systematic literature review.” Computer Methods and
Programs in Biomedicine, 240: 107681.

Pržulj, N., Corneil, D., and Jurisica, I. (2004). ”Modeling interactome:
scale-free or geometric?” Bioinformatics, 20: 3508-3515.

Pudlo, P., Marin, J.-M., Estoup, A., Cornuet, J.-M., Gautier, M., and Robert,
C. (2016). ”Reliable ABC model choice via random forests.” Bioinformatics, 32:
859-866.

Qu, L. and Lu, Y. (2021). ”Copula density estimation by finite mixture of
parametric copula densities.” Communications in Statistics, Simulation and
Computation, 50: 3315-3337.

Raval, R. and Ray, A. (2013). Introduction to Biological Networks. Andover,
U.K.: Taylor & Francis.

Raynal, L., Chen, S., Mira, A., and Onnela, J. (2022). ”Scalable Approximate
Bayesian Computation for growing network models via extrapolated and sampled
summaries.” Bayesian Analysis, 17: 165-192.

Raynal, L., Marin, J., Pudlo, P., Ribatet, M., Robert, C., and Estoup, A.
(2018). ”ABC random forests for Bayesian parameter inference.” Bioinformatics,
35: 1720-1728.

Rosenblatt, J. and Nadler, B. (2016). ”On the optimality of averaging in
distributed statistical learning.” Information and Inference: A Journal of the
IMA, 5: 379-404.

Rossi, R. and Ahmed, N. (2015). ”The network data repository with interactive
graph analytics and visualization.” In Proceedings of the Twenty-Ninth AAAI
Conference on Artificial Intelligence.

Rubio, F. and Johansen, A. (2013). ”A simple approach to maximum intractable
likelihood estimation.” Electronic Journal of Statistics, 7: 1632-1654.

Schmid, C. and Desmarais, B. (2017). ”Exponential random graph models with big
networks: Maximum pseudolikelihood estimation and the parametric bootstrap.”
In 2017 IEEE International Conference on Big Data, 116-121.

Schmid, C. and Hunter, D. (2023). ”Computing pseudolikelihood estimators for
exponential-family random graph models.” Journal of Data Science, 21: 295-309.

Schweinberger, M., Krivitsky, P., Butts, C., and Stewart, J. (2020).
”Exponential family models of random graphs: Inference infinite, super and
infinite population scenarios.” Statistical Science, 35: 627-662.

Sisson, S., Fan, Y., and Beaumont, M. (2018). Handbook of Approximate Bayesian
Computation. Boca Raton, FL: Chapman and Hall/CRC Press.

Sklar, M. (1959). ”Fonctions de repartition an dimensions et leurs marges.”
Public Institute Statistics University Paris, 8: 229-231.

Smith, M. and Klein, N. (2021). ”Bayesian inference for regression copulas.”
Journal of Business and Economic Statistics, 39: 712-728.

Smith, M. and Vahey, S. (2016). ”Asymmetric forecast densities for U.S.
macroeconomic variables from a Gaussian copula model of cross-sectional and
serial dependence.” Journal of Business and Economic Statistics, 34: 416-434.

Snijders, T. (2002). ”Markov chain Monte Carlo estimation of exponential
random graph models.” Journal of Social Structure, 3: 1-40.

Snijders, T. (2011). ”Statistical models for social networks.” Annual Review
of Sociology, 37: 131-153.

Snijders, T., Pattison, P., Robins, G., and Handcock, M. (2006). ”New
specifications for exponential random graph models.” Sociological Methodology,
36: 99-153.

Solé, R., Pastor-Satorras, R., Smith, E., and Kepler, T. (2002). ”A model of
large-scale proteome evolution.” Advances in Complex Systems, 5: 43-54.

Song, P. (2000). ”Multivariate dispersion models generated from Gaussian
copula.” Scandinavian Journal of Statistics, 27: 305-320.

Stewart, J., Schweinberger, M., Bojanowski, M., and Morris, M. (2019).
”Multilevel network data facilitate statistical inference for curved ERGMs
with geometrically weighted terms.” Social Networks, 59: 98-119.

Strauss, D. and Ikeda, M. (1990). ”Pseudolikelihood estimation for social
networks.” Journal of the American Statistical Association, 85: 204-212.

Sunnåker, M., Busetto, A., Numminen, E., Corander, J., Foll, M., and Dessimoz,
C. (2013). ”Approximate Bayesian Computation.” PLOS Computational Biology, 9:
1-10.

Tantardini, M., Ieva, F., Tajoli, L., and Piccardi, C. (2019). ”Comparing
methods for comparing networks.” Scientific Reports, 9: 1-19.

Thiemichen, S., Friel, N., Caimo, A., and Kauermann, G. (2016). ”Bayesian
exponential random graph models with nodal random effects.” Social Networks,
46: 11-28.

Varin, C., Reid, N., and Firth, D. (2011). ”An overview of composite
likelihood methods.” Statistica Sinica, 21: 5-42.

Vazquez, A. (2003). ”Growing network with local rules: Preferential
attachment, clustering hierarchy, and degree correlations.” Physical Review E,
67: 056104.

Vazquez, A., Flammini, A., Maritan, A., and Vespignani, A. (2003). ”Modeling
of protein interaction networks.” ComPlexUs, 1: 38-44.

Wasserman, S. and Faust, K. (1994). Social Network Analysis: Methods and
Applications. U.K.: Cambridge University Press.

Wasserman, S. and Pattison, P. (1996). ”Logit models and logistic regressions
for social networks: I. An introduction to Markov graphs and p*.”
Psychometrika, 61: 401-425.

Watts, D. and Strogatz, S. (1998). ”Collective dynamics of small-world
networks.” Nature, 393: 440-442.

Wei, Z., Kim, S., and Kim, D. (2016). ”Multivariate Skew Normal Copula for
Nonexchangeable Dependence.” Procedia Computer Science, 91: 141-150.

Yang, J. and Leskovec, J. (2012). ”Defining and evaluating network communities
based on ground-truth.” In Proceedings of the ACM SIGKDD Workshop on Mining
Data Semantics, MDS ’12. New York, NY, USA: Association for Computing
Machinery.

Yaveroğlu, Ö., Malod-Dognin, N., Davis, D., Levnajic, Z., Janjic, V.,
Karapandza, R., Stojmirovic, A., and Pržulj, N. (2014). ”Revealing the hidden
language of complex networks.” Scientific reports, 4: 4547.

Yaveroğlu, Ö., Milenković, T., and Pržulj, N. (2015). ”Proper evaluation of
alignment free network comparison methods.” Bioinformatics, 31: 2697-2704.

Yildirim, S., Singh, S., Dean, T., and Jasra, A. (2015). ”Parameter estimation
in hidden Markov Models with intractable likelihoods using sequential Monte
Carlo.” Journal of Computational and Graphical Statistics, 24: 846-865.

Yoshiba, T. (2018). ”Maximum likelihood estimation of skew-t copulas with its
applications to stock returns.” Journal of Statistical Computation and
Simulation, 88: 2489-2506.

Zhu, W., Marin, J., and Leisen, F. (2016). ”A Bootstrap likelihood approach to
Bayesian computation.” Australian and New Zealand Journal of Statistics, 58:
227-244.