Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
40,102 characters · 12 sections · 29 citation commands
Deterministic, quenched and annealed parameter estimation for heterogeneous network models
\email{[email removed]}
\pacs{89.75.Fb; 02.50.Tt; 89.65.Gh}
Over the last twenty years, the growth of network science has impacted several disciplines by establishing new, empirical facts about the structural properties of the related systems. Prominent examples are provided by economics and finance: the growing availability of data has motivated researchers to explore and model the architecture of cryptocurrencies Vallarano2020, interbank networks Bardoscia2021, production networks Ialongo2022 and trading networks Garlaschelli2004,Schweitzer2009,Fronczak2012b,Herman2022.
Modelling the establishment of a connection and the corresponding weight simultaneously poses a serious challenge. Econometrics prescribes to estimate binary and weighted parameters either separately, within the context of hurdle models Mullahy, or jointly, within the context of zero-inflated models Burger2009; in both cases, the Gravity Model specification Tinbergen1962 $\langle w_{ij}\rangle_\text{GM}=f(\omega_i,\omega_j,d_{ij}|\underline{\phi})=e^{\rho}(\omega_i\omega_j)^{\alpha}d_{ij}^{\gamma}$ - where $\omega_i\equiv{\text{GDP}_i}/{\overline{\text{GDP}}}$ is the GDP of country $i$ divided by the arithmetic mean of the GDPs of all countries, $d_{ij}$ is the geographic distance between the capitals of countries $i$ and $j$ and $\underline{\phi}\equiv(\rho,\alpha,\gamma)$ is the vector of parameters defining the Gravity Model specification - is interpreted as the expected value of a probability distribution whose functional form is arbitrary. On the other hand, the approach rooted in statistical physics constructs maximum-entropy distributions, constrained to satisfy certain network properties Shannon,Jaynes1957a,Cover,DSBook,Cimini2019.
In a couple of recent, companion papers Marzio2022,Marzio2023 the two, aforementioned approaches have been integrated within the framework induced by the constrained optimisation of the Kullback-Leibler (KL) divergence Kullback1951. In particular, two, broad classes of models have been constructed, i.e. the integrated and conditional ones, defined by different, probabilistic rules to place links, load them with weights and turn them into properly econometric prescriptions. For what concerns integrated models, the first, two rules follow from a single, constrained optimisation of the KL divergence Garlaschelli2009; for what concerns conditional models, the two rules are disentangled and the functional form of the weight distribution follows from a conditional, optimisation procedure Parisi2020. Still, the prescriptions adopted by the two approaches to carry out the estimation of the parameters entering into the definition of each model differ.
The present contribution is devoted to comparing these recipes in the case of continuous, conditional network models defined by both homogeneous and heterogeneous constraints.
The functional form of continuous, conditional network models can be identified through the constrained minimisation of the KL divergence of a distribution $Q$ from a prior distribution $R$, i.e.
where $\mathbf{W}$ is one of the possible values of a continuous random variable, $\mathbb{W}$ is the set of possible values that $\mathbf{W}$ can take, $Q(\mathbf{W})$ is the (multivariate) probability density function to be estimated and $R(\mathbf{W})$ plays the role of prior distribution, whose divergence from $Q(\mathbf{W})$ must be minimised: in our setting, $\mathbf{W}$ represents an entire network whose weights, now, obey the property $w_{ij}\in\mathbb{R}^+_0$, $\forall\:i<j$. Such an optimisation scheme embodies the so-called Minimum Discrimination Information Principle Marzio2022,Marzio2023, implementing the idea that, as new information becomes available, an updated distribution $Q(\mathbf{W})$ should be chosen in order to make its discrimination from the prior distribution $R(\mathbf{W})$ as hard as possible.
Let us, now, separate both the prior and the posterior distribution into a purely binary part and a conditional, weighted one; the positions $Q(\mathbf{W})=P(\mathbf{A})Q(\mathbf{W}|\mathbf{A})$ and $R(\mathbf{W})=T(\mathbf{A})R(\mathbf{W}|\mathbf{A})$, where $\mathbf{A}$ denotes the binary projection of the weighted network $\mathbf{W}$ (i.e. $\Theta[\mathbf{W}]=\mathbf{A}$), $T(\mathbf{A})$ represents the binary prior and $R(\mathbf{W}|\mathbf{A})$ represents the conditional, weighted prior, lead the KL divergence to be re-writable as
i.e. as a sum of the two addenda
In what follows, we will deal with completely uninformative priors, a choice that amounts at considering the (somehow, simplified) expression
i.e. `minus' the joint entropy, where
is the Shannon entropy of the probability distribution describing the binary projection of the network structure DSBook,Cimini2019 and
is the conditional Shannon entropy of the probability distribution describing the weighted network structure Marzio2022,Marzio2023,Parisi2020. Notice that, when continuous models are considered, $S(\overline{Q}|P)$ is defined by a sum running over all the binary configurations within the ensemble $\mathbb{A}$ and an integral over all the weighted configurations that are compatible with each, specific, binary structure, i.e. $\mathbb{W}_\mathbf{A}=\{\mathbf{W}:\Theta[\mathbf{W}]=\mathbf{A}\}$. For a more detailed discussion, see Appendix A.
The functional form of $P(\mathbf{A})$ can be determined by carrying out the usual, constrained maximisation of Shannon entropy DSBook,Cimini2019; remarkably, any set of (binary) constraints considered in the present paper will lead to the same expression for $P(\mathbf{A})$, i.e. $P(\mathbf{A})=\prod_{i<j}p_{ij}^{a_{ij}}(1-p_{ij})^{1-a_{ij}}$ with $p_{ij}=x_{ij}/(1+x_{ij})$: specifically, the position $x_{ij}\equiv x$ individuates the Undirected Binary Random Graph Model (UBRGM), the position $x_{ij}\equiv x_ix_j$ individuates the Undirected Binary Configuration Model (UBCM) and the position $x_{ij}\equiv \delta\omega_i\omega_j$ individuates the Logit Model (LM) Caldarelli2003.
On the other hand, the functional form of $Q(\mathbf{W}|\mathbf{A})$ can be determined by carrying out the constrained maximisation of $S(\overline{Q}|P)$, the set of constraints being, now,
while the first condition ensures the normalisation of the probability distribution, the vector $\{C_\alpha(\mathbf{W})\}$ represents the proper set of weighted constraints. The distribution induced by such an optimisation problem reads
if $\mathbf{W}\in\mathbb{W}_\mathbf{A}$ and $0$ otherwise. While the Hamiltonian $H(\mathbf{W})=\sum_\alpha\psi_\alpha C_\alpha(\mathbf{W})$ lists the constraints, the quantity at the denominator is the partition function, conditional on the fixed topology $\mathbf{A}$ Parisi2020.
For mathematical convenience, in what follows we will consider separable Hamiltonians, i.e. functions that can be written as sums of node pair-specific Hamiltonians: $H(\mathbf{W})=\sum_{i<j}H_{ij}(w_{ij})$; this choice leads to the result
(with $m_{ij}$ being the pair-specific, minimum weight allowed by a given model and $\zeta_{ij}$ being the corresponding partition function), irrespectively from the specific, functional form of $H_{ij}(w_{ij})$ Marzio2023. For a more detailed discussion, see Appendix B.
Several, alternative recipes are viable to estimate the parameters entering into the definition of continuous, conditional network models.
The simplest one prescribes to consider the traditional likelihood function
with $\mathbf{W}^*$ ($\mathbf{A}^*$) being the empirical, weighted (binary) adjacency matrix; its maximisation allows the parameters entering into the definition of the purely topological distribution and those entering into the definition of the conditional, weighted one to be estimated in a totally disentangled fashion Marzio2023. In fact, maximising
with respect to the unknown parameters leads us to find the vector of values $\underline{\psi}^*$ satisfying the vector of relationships
which stands for the set of relationships $\langle C_\alpha\rangle_{\mathbf{A}^*}(\underline{\psi}^*)\equiv C_\alpha^*$, $\forall\:\alpha$, each one equating the model-induced, average value of the corresponding constraint to its empirical value, marked with an asterisk.
This first approach to parameter estimation can be named as `deterministic', to stress that $\mathbf{A}^*$ is considered as not being subject to variation; otherwise stated, this recipe - which is the most common in econometrics - prescribes to estimate the parameters entering into the definition of the conditional, weighted probability distribution by assuming the network topology to be fixed.
Topology, however, is a random variable itself, obeying the probability distribution $P(\mathbf{A})$. As a consequence, the `deterministic' recipe for parameter estimation could lead to inconsistencies, should the description of $\mathbf{A}^*$ provided by $P(\mathbf{A})$ be not accurate. The variability induced by $P(\mathbf{A})$ can be properly accounted for by considering the generalised likelihood Parisi2020
whose maximisation leads us to find the vector of values $\underline{\psi}^*$ satisfying the vector of relationships
which stands for the set of relationships $\langle C_\alpha\rangle(\underline{\psi}^*)\equiv C_\alpha^*$, $\forall\:\alpha$. Taking this average is conceptually similar to taking the `annealed' average in physics: parameter estimation is carried out while random variables - again, the entries of the adjacency matrix - are left to vary.
Interestingly, the `deterministic' recipe is a special case of the `annealed' recipe since the former can be recovered by posing $P(\mathbf{A})\equiv\delta_{\mathbf{A},\mathbf{A}^*}$: in this case, in fact,
similarly, $\sum_{\mathbf{A}\in\mathbb{A}}\delta_{\mathbf{A},\mathbf{A}^*}\langle\mathbf{C}\rangle_{\mathbf{A}}(\underline{\psi}^*)=\langle\mathbf{C}\rangle_{\mathbf{A}^*}(\underline{\psi}^*)=\mathbf{C^*}$.
A viable alternative to properly account for the variability induced by $P(\mathbf{A})$ is that of reversing the two operations of `likelihood maximisation' and `ensemble averaging': in other words, one can 1) numerically sample the ensemble of configurations induced by $P(\mathbf{A})$, 2) maximise the likelihood $\ln Q(\mathbf{W}^*|\mathbf{A})$ for each, generated network, 3) take the average of the resulting set of parameters, according to the formula
the estimation of the $\alpha$-th parameter being assumed to coincide with the average $\langle\psi_\alpha^*\rangle$.
Taking this average is conceptually similar to taking the `quenched' average in physics: random variables - in the specific case, the entries of the adjacency matrix - are frozen, parameter estimation is carried out and, only at the end, the values of the parameters are averaged over the ensemble of configurations induced by $P(\mathbf{A})$.\\
As our models inherit their functional form from the constrained minimisation of the KL divergence, each parameter controls for a specific constraint: when employing the `deterministic' recipe, such a circumstance makes each parameter configuration-dependent; when employing either the `annealed' or the `quenched' recipe, instead, accounting for the variability of a network structure induces a sort of `loss of memory' about its empirical, purely topological details.
In order to test if the `deterministic', `annealed' and `quenched' prescriptions lead to the same estimation, let us focus on a number of variants of the Conditional Exponential Model (CEM), induced by the positions $H_{ij}^\text{CEM}=\beta_{ij}w_{ij}$ and $\zeta_{ij}^\text{CEM}=\beta_{ij}^{-1}$:
naturally, $q_{ij}(w_{ij}=0|a_{ij}=0)=1$ (i.e. if nodes $i$ and $j$ are not connected, the weight of the corresponding link is zero with probability equal to one) and $q_{ij}(w_{ij}>0|a_{ij}=1)=\beta_{ij}e^{-\beta_{ij}w_{ij}}$.
In what follows, we will consider three, different instances of $p_{ij}=x_{ij}/(1+x_{ij})$, corresponding to
Let us start by considering the `scalar' or homogeneous variant of the CEM, defined by the position $\beta_{ij}\equiv\beta$, $\forall\:i<j$.
In this case, the `deterministic' recipe for parameter estimation prescribes to maximise the likelihood
where $W(\mathbf{W}^*)\equiv W^*=\sum_{i<j}w_{ij}^*$ and whose optimisation leads to the expression $\beta=L^*/W^*$. The `annealed' recipe prescribes to maximise the likelihood
whose optimisation leads to the expression $\beta=\langle L\rangle/W^*$. The `quenched' recipe, on the other hand, prescribes to calculate the average
since, now, $\beta(\mathbf{A})=L(\mathbf{A})/W^*$.
In the case of the `scalar' variant of the CEM, the estimations coincide for any null model preserving the total number of links, i.e. ensuring that $\langle L\rangle=L^*$, regardless of the network density. Such a result is confirmed by Fig. (ref) where each recipe has been implemented on the WTW, by adopting the distributions induced by the UBRGM (blue), the UBCM (green) and the LM (red). Specifically, the `deterministic' estimation (black, solid line) and the `annealed' estimations (blue, green and red, solid lines) overlap; moreover, each `annealed' estimation overlaps with the the corresponding, `quenched' estimation, i.e. the average value of the related, `quenched' distribution (blue, green and red, dash-dotted lines).
In the case of the UBRGM-induced, homogeneous version of the CEM, the `quenched' distribution of the parameter $\beta(\mathbf{A})=L(\mathbf{A})/W^*$ `inherits' the distribution of the total number of links, i.e. $L\sim\text{Bin}(N(N-1)/2,p)$, with $p=2L^*/N(N-1)$: more precisely, $W\beta\sim\text{Bin}(N(N-1)/2,p)$; analogously for the UBCM- and the LM-induced, homogeneous versions of the CEM - the only difference being that, now, $L$ obeys two, different, Poisson-Binomial distributions.
Let us, now, consider the `vector' or weakly heterogeneous variant of the CEM, defined by the position $\beta_{ij}\equiv\beta_i+\beta_j$, $\forall\:i<j$.
In this case, the `deterministic' recipe for parameter estimation prescribes to maximise the likelihood
where $s_i(\mathbf{W}^*)\equiv s_i^*=\sum_{j(\neq i)}w_{ij}^*$ and whose optimisation requires to solve the system of equations
The `annealed' recipe, instead, prescribes to maximise the likelihood
whose optimisation requires to solve the system of equations
(notice that both the `deterministic' and the `annealed' version of the `vector' variant of the CEM are alternative instances of the so-called $\text{CReM}_\text{A}$, introduced in Parisi2020). The `quenched' recipe, on the other hand, requires to solve the system of equations $\langle\beta_i\rangle=\sum_{\mathbf{A}\in\mathbb{A}}P(\mathbf{A})\beta_i(\mathbf{A})$, $\forall\:i$ which no longer have an explicit expression. Devising some sort of approximation is, however, possible. Let us start by re-writing eq. (ref) as
and consider the node whose coefficient is the largest one. This allows us to write $\beta_i\simeq\sum_{j(\neq i)}p_{ij}/s_i^*=\langle k_i\rangle/s_i^*$: in case we implemented the UBRGM, we would obtain $\beta_i(\mathbf{A})\simeq 2L(\mathbf{A})/Ns_i^*$, hence expecting the `quenched' distribution of $Ns_i^*\beta_i/2$ to coincide with $\text{Bin}(N(N-1)/2,p)$; if, on the other hand, we implemented the UBCM, we would obtain $\beta_i(\mathbf{A})\propto k_i(\mathbf{A})/s_i^*$, hence expecting the `quenched' distribution of $s_i^*\beta_i$ to obey a Poisson-Binomial. Again, the estimations coincide for any null model preserving the structural properties characterising the binary recipe implemented.
More generally, the mutual relationships between the estimations provided by the three recipes are node-dependent (see Fig. (ref), illustrating the case-study of node 166 of the WTW and Fig. (ref) in Appendix C): in general, however, each `annealed' estimation overlaps with the average value of the related `quenched' distribution. Moreover, the `deterministic' estimation is very close to the UBCM-induced, `annealed' one; such a result is a consequence of the accurate description of the empirical network topology provided by the UBCM - in fact, much more accurate than the ones provided by the UBRGM and the LM: indeed, the better the approximation $p_{ij}\simeq a_{ij}$, $\forall\:i<j$, the closer the `annealed' estimation to the `deterministic' one.
This is even more evident when considering the `tensor' variant of the CEM, in which case the three optimisation procedures lead to the expressions $\beta_\text{det}=a_{ij}^*/\hat{w}_{ij}$, $\forall\:i<j$ and $\beta_\text{ann}=\langle\beta\rangle_\text{que}=p_{ij}/\hat{w}_{ij}$, $\forall\:i<j$ - with $\hat{w}_{ij}$ representing an estimate of the empirical weight $w_{ij}^*$; if, however, $\hat{w}_{ij}\equiv w_{ij}^*$, $\forall\:i<j$ then, for consistency, $p_{ij}\equiv a_{ij}^*$ and the three recipes coincide.
As a third case-study, let us focus on the `econometric' variant of the CEM, defined by posing $\beta_{ij}\equiv\beta_0+z_{ij}^{-1}$, $\forall\:i<j$, where $z_{ij}\equiv e^{\rho}(\omega_i\omega_j)^\alpha d_{ij}^\gamma$ represents the Gravity Model specification traditionally employed to analyse undirected, weighted, trade networks and $\beta_0$ is a structural parameter to be tuned in order to ensure that $\langle W\rangle=W^*$. In this case, the `deterministic' recipe for parameter estimation prescribes to maximise the likelihood
whose optimisation requires to solve the system of equations
The `annealed' recipe, instead, prescribes to maximise the likelihood
whose optimisation requires to solve the system of equations
The `quenched' recipe, on the other hand, requires to solve the system of equations $\langle\beta_{0}\rangle=\sum_{\mathbf{A}\in\mathbb{A}}P(\mathbf{A})\beta_{0}(\mathbf{A})$ and $\langle\underline{\phi}\rangle=\sum_{\mathbf{A}\in\mathbb{A}}P(\mathbf{A})\underline{\phi}(\mathbf{A})$ which no longer have an explicit expression.
Figures (ref) and (ref) in Appendix C illustrate the case-study of the WTW: although the `quenched' distributions induced by the three, binary recipes are characterised by different shapes that may overlap (as in the case of the parameters $\rho$ - under the UBRGM-induced and UBCM-induced binary recipes - and $\gamma$ - under all, binary recipes) or not (as in the case of the parameters $\beta_0$ and $\alpha$), `annealed' and `quenched' estimations always coincide (the only, small discrepancy being observable for the parameter $\beta_0$, under the UBRGM-induced, binary recipe). The `deterministic' estimation, instead, is compatible with the other, two ones only for the parameter $\alpha$, under the UBCM-induced, binary recipe.\\
Sparse networks deserve a separate discussion. The results concerning the homogeneous and econometric variant of the BLN, defined by posing $\beta_{ij}\equiv\beta_0+z_{ij}^{-1}$, $\forall\:i<j$, with $z_{ij}\equiv e^{\rho}(s_i s_j)^\alpha$, are analogous to the ones shown for the WTW - in the latter case, the `annealed' estimates of $\beta_0$, $\rho$ and $\alpha$ are very close to their `quenched' counterparts, the relative error $\text{RE}=|(\phi_i^\text{ann}-\phi_i^\text{que})/\phi_i^\text{ann}|$ amounting at $\simeq 10^{-3}$ for $\beta_0$ and $\simeq 10^{-4}$ for $\rho$, $\alpha$. On the contrary, these conclusions no longer hold true when the weakly heterogeneous variant of the CEM is considered: in this case, in fact, carrying out the `quenched' approach can lead to binary configurations with disconnected nodes, a circumstance that impairs the correct estimation of the corresponding parameters; carrying out the `annealed' estimation, instead, remains a feasible task.
The present contribution focuses on three recipes for estimating the parameters entering into the definition of statistical network models, i.e. the `deterministic', `annealed' and `quenched' ones. In order to implement them, we have considered several variants of the CEM, i.e. the homogeneous one (defined by one, global parameter), the weakly heterogeneous one (defined by $N$, local parameters) and the econometric one (defined by four, global parameters), each one combined with three, different recipes for estimating the network topology (i.e. the UBRGM, the UBCM and the LM).
The `deterministic' recipe, routinely employed in econometrics to determine the so-called hurdle models Mullahy, prescribes to estimate the parameters associated to the weighted constraints on the empirical realisation of the network topology. Since it considers $\mathbf{A}^*$ as not being subject to variation, its use is recommended whenever $\text{Var}[a_{ij}]=p_{ij}(1-p_{ij})\simeq 0$ or, equivalently, $p_{ij}\simeq a_{ij}$, $\forall\:i<j$, i.e. whenever the binary random variables can be safely considered as deterministic or, more in general, whenever their (scale of) variation is negligible with respect to the (scale of) variation of the weighted random variables.
Accounting for such a variability in a fully consistent manner can be achieved upon adopting either the `annealed' recipe (according to which parameters are estimated on the average network topology) or the `quenched' recipe (according to which parameters are, first, estimated on a large number of binary configurations and, then, averaged); the main difference between these procedures lies in the order in which the two operations of `averaging' (of the entries of the binary adjacency matrix) and `maximisation' (of the related likelihood function) are taken. Interestingly, no variant of the CEM is sensitive to this choice (neither the purely structural ones nor the `econometric' one); while, however, the coincidence of the `annealed' and `quenched' estimates for purely structural models can be explicitly verified, this is no longer true when the `econometric' variant is considered: in this case, in fact, one can proceed only numerically.
This evidence reveals the main limitation of the `quenched' approach, i.e. the need of resorting upon an explicit sampling of the chosen, binary ensemble. As any `good' sampling algorithm must lead to a faithful representation of the parent distribution, we are left with the following question: is this always guaranteed, in all cases of interest to us?
This seems to be the case for dense networks. As shown in Tiziano2015, a study of the coefficient of variation of the constraints defining the `vector' variant of the CEM (i.e. the ratio between standard deviation and expected value of each degree) reveals it to vanish in the asymptotic limit: in other words, the fluctuations affecting each degree vanish, a result guaranteeing that the degree sequence of any configuration in the ensemble remains `close enough' to the empirical one.
When sparse networks are, instead, considered, the coefficient of variation of the constraints defining the `vector' variant of the CEM remains finite in the asymptotic limit: in other words, the fluctuations affecting each degree do not vanish, a result implying that the degree sequence of any configuration in the ensemble may largely differ from the empirical one; to provide a concrete example, nodes whose empirical degree is `small' may disconnect, hence inducing the resolution of a system of equations which is not even compatible with the set of constraints defining the original problem. Overcoming such a limitation implies quantifying the bias affecting the estimates in cases like these: although possible, calculations of this kind are far beyond the scope of the present paper.\\
Overall, then, two alternatives exist to overcome the main limitation of the `deterministic' estimation recipe, i.e. that of ignoring the variety of structures that are compatible with a given probability distribution $P(\mathbf{A})$, namely the `annealed' and `quenched' ones. As the `quenched' recipe requires an explicit sampling the ensemble - potentially leading to inconsistent estimates for sparse configurations - we believe the `annealed' one to represent the better alternative, 1) being unbiased by definition, 2) being convenient from a numerical point of view, 3) reducing to the `deterministic' recipe in case the empirical configuration is not subject to variation.
SoBigData.it receives funding from European Union – NextGenerationEU – National Recovery and Resilience Plan (Piano Nazionale di Ripresa e Resilienza, PNRR) – Project: “SoBigData.it – Strengthening the Italian RI for Social Mining and Big Data Analytics” – Prot. IR0000013 – Avviso n. 3264 del 28/12/2021. This work is also supported by PNRR-M4C2-Investimento 1.3, Partenariato Esteso PE00000013 - `FAIR-Future Artificial Intelligence Research' - Spoke 1 `Human-centered AI', funded by the European Commission under the NextGeneration EU programme and by the project `Network analysis of economic and financial resilience', Italian DM n. 289, 25-03-2021 (PRO3 Scuole) CUP D67G22000130001. DG acknowledges support from the Dutch Econophysics Foundation (Stichting Econophysics, Leiden, the Netherlands) and the Netherlands Organization for Scientific Research (NWO/OCW). MDV acknowledges support from the European Union ERC-2018-ADG Grant Agreement n. 834756, `XAI: Science and technology for the explanation of AI decision making'. MDV and DG also acknowledge support from the `Programma di Attività Integrata' (PAI) project `Prosociality, Cognition and Peer Effects' (Pro.Co.P.E.), funded by IMT School for Advanced Studies Lucca.