Bayesian model selection for exponential random graph models via adjusted pseudolikelihoods

TL;DR

This paper introduces an adjusted pseudolikelihood approach for Bayesian model selection in exponential random graph models, enabling more efficient and reliable inference in complex network analysis.

Contribution

It proposes a novel adjustment method for pseudolikelihoods to approximate likelihoods accurately, facilitating Bayesian model comparison in doubly-intractable models.

Findings

The adjusted pseudolikelihood provides evidence estimates comparable to existing methods.

The proposed method reduces computational cost significantly.

Empirical results demonstrate improved inference quality in network models.

Abstract

Models with intractable likelihood functions arise in areas including network analysis and spatial statistics, especially those involving Gibbs random fields. Posterior parameter es timation in these settings is termed a doubly-intractable problem because both the likelihood function and the posterior distribution are intractable. The comparison of Bayesian models is often based on the statistical evidence, the integral of the un-normalised posterior distribution over the model parameters which is rarely available in closed form. For doubly-intractable models, estimating the evidence adds another layer of difficulty. Consequently, the selection of the model that best describes an observed network among a collection of exponential random graph models for network analysis is a daunting task. Pseudolikelihoods offer a tractable approximation to the likelihood but should be treated with caution because they can lead to an unreasonable inference. This paper specifies a method to adjust pseudolikelihoods in order to obtain a reasonable, yet tractable, approximation to the likelihood. This allows implementation of widely used computational methods for evidence estimation and pursuit of Bayesian model selection of exponential random graph models for the analysis of social networks. Empirical comparisons to existing methods show that our procedure yields similar evidence estimates, but at a lower computational cost.

Tables

Table 1. Table 1: Zachary karate club - CPU time in minutes for each adjustment phase.

Model	Mode	Curvature	Magnitude	Total
${ℳ}_1$	0.336	0.002	2.688	3.026
${ℳ}_3$	0.417	0.002	2.874	3.293

Table 2. Table 2: Zachary karate club - Average and standard deviation values of log-marginal likelihood, Bayes factor estimates and CPU time in minutes from thirty independent experiments, based on the unadjusted and the fully adjusted pseudolikelihood function. The CPU time corresponds to the total computational time required to apply an algorithm to both models of interest.

Adjustment	Method	$\log \pi \left(y\mid {ℳ}_1\right)$	$\log \pi \left(y\mid {ℳ}_3\right)$	$B{F}_{13}$	CPU
(a) No	Chib & Jeliazkov	-217.197 (0.01)	-219.842 (0.01)	14.088 (0.04)	14.20
	Stepping stones	-216.805 (0.23)	-219.520 (0.26)	15.809 (4.71)	46.17
	Power posteriors - TI	-216.798 (0.24)	-219.523 (0.27)	16.041 (5.01)	46.17
	Power posteriors - CTI	-217.077 (0.05)	-219.726 (0.04)	14.173 (0.81)	46.17
(b) Yes	Chib & Jeliazkov	-219.007 (0.01)	-221.766 (0.01)	15.776 (0.06)	20.79
	Stepping stones	-218.765 (0.15)	-221.524 (0.20)	16.192 (3.79)	57.90
	Power posteriors - TI	-218.763 (0.16)	-221.525 (0.21)	16.302 (4.00)	57.90
	Power posteriors - CTI	-218.967 (0.02)	-221.716 (0.02)	15.631 (0.42)	57.90
	Population exchange	-218.954 (0.02)	-221.703 (0.02)	15.629 (0.35)	1749.60
	Auto-RJ	-	-	13.177 (0.24)	4777.86

Table 3. Table 3: Teenage Friends and Lifestyle - CPU time in minutes for each adjustment phase.

Model	Mode	Curvature	Magnitude	Total
${ℳ}_1$	0.291	0.002	5.204	5.497
${ℳ}_2$	0.297	0.002	5.094	5.393

Table 4. Table 4: Teenage Friends and Lifestyle - Average and standard deviation values of log-marginal likelihood, Bayes factor estimates and CPU time in minutes from thirty independent experiments, based on the unadjusted and the fully adjusted pseudolikelihood function. The CPU time corresponds to the total computational time required to apply an algorithm to both models of interest.

Adjustment	Method	$\log \pi \left(y\mid {ℳ}_1\right)$	$\log \pi \left(y\mid {ℳ}_2\right)$	$B{F}_{21}$	CPU
(a) No	Chib & Jeliazkov	-152.010 (0.01)	-153.153 (0.01)	0.319 (0.01)	18.75
	Stepping stones	-152.204 (0.14)	-153.316 (0.20)	0.340 (0.09)	52.41
	Power posteriors - TI	-152.206 (0.15)	-153.312 (0.21)	0.343 (0.09)	52.41
	Power posteriors - CTI	-152.215 (0.01)	-153.362 (0.01)	0.318 (0.01)	52.41
(b) Yes	Chib & Jeliazkov	-233.740 (0.01)	-233.621 (0.01)	1.127 (0.01)	28.13
	Stepping stones	-233.965 (0.12)	-233.801 (0.20)	1.212 (0.28)	67.44
	Power posteriors - TI	-233.967 (0.12)	-233.797 (0.20)	1.221 (0.29)	67.44
	Power posteriors - CTI	-233.975 (0.01)	-233.847 (0.01)	1.137 (0.02)	67.44
	Population exchange	-233.934 (0.01)	-233.807 (0.01)	1.135 (0.02)	7026.06
	Auto-RJ	-	-	1.186 (0.02)	2977.08

Equations

$$f(y\mid \theta )=\frac{q(y\mid \theta )}{z(\theta )}$$

$$\pi (\theta \mid y)=\frac{f(y\mid \theta )p(\theta )}{\pi (y)}=\frac{f(y\mid \theta )p(\theta )}{{\int }_{\Theta }f(y\mid \theta )p(\theta )d\theta }$$

$$\pi ({\theta }_m,{\mathcal{M}}_m\mid y)\propto f(y\mid {\theta }_m,{\mathcal{M}}_m)p({\theta }_m\mid {\mathcal{M}}_m)p({\mathcal{M}}_m).$$

$$\pi (y\mid {\mathcal{M}}_m)={\int }_{{\Theta }_m}f(y\mid {\theta }_m,{\mathcal{M}}_m)p({\theta }_m\mid {\mathcal{M}}_m)\mathrm{d}{\theta }_m,$$

$$\pi ({\mathcal{M}}_m\mid y)=\frac{\pi (y\mid {\mathcal{M}}_m)p({\mathcal{M}}_m)}{{\sum }_{j\in M}\pi (y\mid {\mathcal{M}}_j)p({\mathcal{M}}_j)}$$

$$B{F}_{m,m^{\prime }}=\frac{\pi (y\mid {\mathcal{M}}_m)}{\pi (y\mid {\mathcal{M}}_{m^{\prime }})}.$$

$$y_{ij}=\{\begin{array}{ll}1,& \text{if }(i,j)\text{ are connected,}\\ 0,& \text{otherwise.}\end{array}$$

$$f(y\mid \theta )=\frac{q(y\mid \theta )}{z(\theta )}=\frac{\exp \left\{{\theta }^{\top }s(y)\right\}}{{\sum }_{y\in \mathcal{Y}}\exp \left\{{\theta }^{\top }s(y)\right\}},{\theta }^{\top }s(y)=\sum _{j=1}^d{\theta }_j{s}_j(y),$$

$$\text{logit}\left\{p(y_{ij}=1\mid y_{-ij},\theta )\right\}={\theta }^{\top }{\delta }_s(y{)}_{ij},$$

$$f_{\text{PL}}(y\mid \theta )=\prod _{\begin{array}{c}i\ne j\\ i<j\end{array}}p(y_{ij}\mid y_{-ij},\theta )=\prod _{\begin{array}{c}i\ne j\\ i<j\end{array}}\frac{p(y_{ij}=1\mid y_{-ij},\theta {)}^{y_{ij}}}{\{1-p(y_{ij}=1\mid y_{-ij},\theta ){\}}^{y_{ij}-1}},$$

$$g:\{\begin{array}{l}\Theta \to \Theta \\ \theta \mapsto {\hat{\theta }}_{MPLE}+W(\theta -{\hat{\theta }}_{MLE}),\end{array}$$

$$\tilde{f}(y\mid \theta )=C\cdot f_{\text{PL}}(y\mid g(\theta )),$$

$${\hat{\theta }}_{MLE}$$

$${\hat{\theta }}_{MPLE}$$

$${\nabla }_{\theta }^2\log f(y\mid \theta ){|}_{{\hat{\theta }}_{MLE}}=W^{\top }{\nabla }_{\theta }^2\log f_{\text{PL}}(y\mid \theta ){|}_{{\hat{\theta }}_{MPLE}}W,$$

$${\nabla }_{\theta }\log f(y\mid \theta )$$

$$=s(y)-\frac{\sum s(y)\text{exp}\left\{{\theta }^{\top }s(y)\right\}}{\sum \text{exp}\left\{{\theta }^{\top }s(y)\right\}}$$

$$=s(y)-{\mathbb{E}}_{y\mid \theta }\left[s(y)\right],$$

$${\nabla }_{\theta }^2\log f(y\mid \theta )$$

$$\displaystyle=-\bigg{\{}\mathbb{E}_{y\mid\theta}\left[s^{2}(y)\right]-\left[\mathbb{E}_{y\mid\theta}\left[s(y)\right]\right]^{2}\bigg{\}}$$

$$=-{\mathbb{V}}_{y\mid \theta }\left[s(y)\right],$$

$$-{\nabla }_{\theta }^2\log f(y\mid \theta ){|}_{{\hat{\theta }}_{MLE}}$$

$$-{\nabla }_{\theta }^2\log f_{\text{PL}}(y\mid \theta ){|}_{{\hat{\theta }}_{MPLE}}$$

$$C=\frac{q(y\mid {\hat{\theta }}_{MLE})\cdot z^{-1}({\hat{\theta }}_{MLE})}{f_{\text{PL}}(y\mid g({\hat{\theta }}_{MLE}))}.$$

$$f(y\mid t_j\theta )=\frac{q(y\mid t_j\theta )}{z(t_j\theta )}=\frac{\exp \{(t_j\theta {)}^{\top }s(y)\}}{{\sum }_{y\in \mathcal{Y}}\exp \{(t_j\theta {)}^{\top }s(y)\}},j\in \{0,\dots ,L\}.$$

$$\frac{z({\hat{\theta }}_{MLE})}{z(0)}=\frac{z(t_L{\hat{\theta }}_{MLE})}{z(t_0{\hat{\theta }}_{MLE})}=\prod _{j=0}^{L-1}\frac{z(t_{j+1}{\hat{\theta }}_{MLE})}{z(t_j{\hat{\theta }}_{MLE})},$$

$$\frac{z(t_{j+1}{\hat{\theta }}_{MLE})}{z(t_j{\hat{\theta }}_{MLE})}={\mathbb{E}}_{y\mid t_j{\hat{\theta }}_{MLE}}\left[\frac{q(y\mid t_{j+1}{\hat{\theta }}_{MLE})}{q(y\mid t_j{\hat{\theta }}_{MLE}))}\right].$$

$$\hat{\frac{z(t_{j+1}{\hat{\theta }}_{MLE})}{z(t_j{\hat{\theta }}_{MLE})}}=\frac{1}{K}\sum _{k=1}^K\frac{q(y_k^{(j)}\mid t_{j+1}{\hat{\theta }}_{MLE})}{q(y_k^{(j)}\mid t_j{\hat{\theta }}_{MLE})}.$$

$$\hat{C}=\frac{q(y\mid {\hat{\theta }}_{MLE})\cdot {\hat{z}}^{-1}({\hat{\theta }}_{MLE})}{f_{\text{PL}}(y\mid g({\hat{\theta }}_{MLE}))},$$

$${\pi }_{\text{PL}}({\theta }_m\mid y,{\mathcal{M}}_m)=\frac{f_{\text{PL}}(y\mid {\theta }_m,{\mathcal{M}}_m)p({\theta }_m\mid {\mathcal{M}}_m)}{{\pi }_{\text{PL}}(y\mid {\mathcal{M}}_m)}=\frac{f_{\text{PL}}(y\mid {\theta }_m,{\mathcal{M}}_m)p({\theta }_m\mid {\mathcal{M}}_m)}{{\int }_{{\Theta }_m}{f}_{\text{PL}}(y\mid {\theta }_m,{\mathcal{M}}_m)p({\theta }_m\mid {\mathcal{M}}_m)\mathrm{d}{\theta }_m}.$$

$$B{F}_{m{m}^{\prime }}^{PL}=\frac{{\pi }_{\text{PL}}(y\mid {\mathcal{M}}_m)}{{\pi }_{\text{PL}}(y\mid {\mathcal{M}}_{m^{\prime }})}=\frac{{\int }_{{\Theta }_m}{f}_{\text{PL}}(y\mid {\theta }_m,{\mathcal{M}}_m)p({\theta }_m\mid {\mathcal{M}}_m)\mathrm{d}{\theta }_m}{{\int }_{{\Theta }_{m^{\prime }}}{f}_{\text{PL}}(y\mid {\theta }_{m^{\prime }},{\mathcal{M}}_{m^{\prime }})p({\theta }_{m^{\prime }}\mid {\mathcal{M}}_{m^{\prime }})\mathrm{d}{\theta }_{m^{\prime }}}.$$

$$\tilde{\pi }(\theta \mid y,{\mathcal{M}}_m)=\frac{\tilde{f}(y\mid {\theta }_m,{\mathcal{M}}_m)p({\theta }_m\mid {\mathcal{M}}_m)}{\tilde{\pi }(y\mid {\mathcal{M}}_m)}=\frac{\tilde{f}(y\mid {\theta }_m,{\mathcal{M}}_m)p({\theta }_m\mid {\mathcal{M}}_m)}{{\int }_{{\Theta }_m}\tilde{f}(y\mid {\theta }_m,{\mathcal{M}}_m)p({\theta }_m\mid {\mathcal{M}}_m)\mathrm{d}{\theta }_m},$$

$${\widetilde{BF}}_{m{m}^{\prime }}=\frac{\tilde{\pi }(y\mid {\mathcal{M}}_m)}{\tilde{\pi }(y\mid {\mathcal{M}}_{m^{\prime }})}=\frac{{\int }_{{\Theta }_m}\tilde{f}(y\mid {\theta }_m,{\mathcal{M}}_m)p({\theta }_m\mid {\mathcal{M}}_m)\mathrm{d}{\theta }_m}{{\int }_{{\Theta }_{m^{\prime }}}\tilde{f}(y\mid {\theta }_{m^{\prime }},{\mathcal{M}}_{m^{\prime }})p({\theta }_{m^{\prime }}\mid {\mathcal{M}}_{m^{\prime }})\mathrm{d}{\theta }_{m^{\prime }}}.$$

Full text

Keywords: Bayes factors, Evidence, Intractable normalising constants.

1 Introduction

Bayesian inference for models that are characterized by an intractable likelihood function has received considerable attention by the statistical community, notably the class of Gibbs random fields (GRFs). Popular examples include the autologistic model (Besag, 1972), used to model the spatial distribution of binary random variables defined on a lattice or grid and the exponential random graph model (ERGM) for social network analysis (Robins et al., 2007). Despite their popularity, posterior parameter estimation for GRFs presents considerable difficulties because the normalising constant $z(\theta)$ of the likelihood density

[TABLE]

is typically intractable for all but trivially small graphs. The posterior distribution defined as

[TABLE]

is termed doubly-intractable because of the intractability of the normalising term of the likelihood model within the posterior and the intractability of the posterior normalising term.

Bayesian model comparison is often achieved by estimating the Bayes factor (Kass and Raftery, 1995), which relies upon the marginal likelihood or model evidence, $\pi(y)$, of each of the competing models. However, for many models of interest with intractable likelihoods such as GRFs, estimation of the marginal likelihood adds another layer of difficulty. This paper addresses this problem in the context of Bayesian model comparison of exponential random graph models.

Related work by Friel (2013) and Everitt et al. (2017) has the same objective as our study, namely to estimate the evidence in the presence of an intractable likelihood normalising constant. Contrary to our method, their proposed algorithms rely heavily on repeated simulations from the likelihood. Friel (2013) devised a "population" version of the exchange algorithm (Møller et al., 2006), however for evidence estimation, it is limited to models with a small number of parameters. Everitt et al. (2017) describe an importance sampling approach for estimating the evidence, which is promising for low-dimensional models. However, when moving to higher dimensional settings their approach makes use of a particle filter to estimate the evidence, which is naturally more computationally demanding.

Motivated by overcoming the intractability of the likelihood in (1), a natural approach is to use composite likelihoods as a plug-in for the true likelihood (Varin et al., 2011). The pseudolikelihood (Besag, 1975) is an antecedent of composite likelihoods and was developed in the context of ERGMs by Strauss and Ikeda (1990). Building on the work of Stoehr and Friel (2015), Bouranis et al. (2017) proposed an alternative approach to Bayesian inference for ERGMs. The replacement of the true likelihood with the pseudolikelihood approximation in Bayes formula yields what is termed a pseudo-posterior distribution, a tractable Bayesian model from which it is straightforward to sample. Bayesian inference based on the pseudolikelihood can be problematic however, as in some cases the posterior mean estimates are biased and the posterior variances are typically underestimated. Bouranis et al. (2017) developed an approach to allow for correction of a sample from the pseudo-posterior distribution so that it is approximately distributed from the target posterior distribution.

While parameter inference based on the adjusted pseudo-posterior distribution yields reasonable results (Bouranis et al., 2017), evidence estimation using this correction procedure is inefficient, a point which is explained in Section 5. Based on this observation, we consider in this paper adjusting the pseudolikelihood directly, as opposed to the pseudo-posterior, and the likelihood function in the model evidence is replaced with this fully adjusted pseudolikelihood. These adjustments involve a correction of the mode, the curvature and the magnitude at the mode of the pseudolikelihood function, as outlined in Figure 1. The crucial point is that this adjusted pseudolikelihood function renders the corresponding posterior distribution amendable to standard evidence estimation methods from the Bayesian toolbox. A non-exhaustive list of such methods includes Chib’s method (Chib, 1995) and its extension (Chib and Jeliazkov, 2001), importance sampling (Liu, 2001), annealed importance sampling (Neal, 2001), bridge sampling (Meng and Wong, 1996) and path sampling/thermodynamic integration (Gelman and Meng, 1998; Lartillot and Phillipe, 2006; Friel and Pettitt, 2008; Calderhead and Girolami, 2009), among others. The tractability of the fully adjusted pseudolikelihood allows for evidence estimation with such methods, thereby allowing Bayesian model selection of exponential random graph models. The method also applies to more general Boltzmann distributions, used in statistical physics.

The outline of the paper is as follows. Section 2 introduces the reader to the concept of Bayesian model comparison. A basic description of exponential random graph models is provided in Section 3. In Section 4 we discuss how to perform the adjustments of the pseudolikelihood for ERGMs with the goal to obtain an approximation of the marginal likelihood and in Section 5 we derive an approximation to the Bayes factor. In Section 6, we assess the efficiency of the marginal likelihood approximation with a Potts model example for spatial analysis (Potts, 1952), where the size of the lattice allows for exact estimation of the marginal likelihood. Detailed ERG model selection experiments are presented in Section 7. We conclude the paper in Section 8 with final remarks and recommendations to practitioners based on accuracy of evidence and Bayes factor estimates and computational speed. The Bergm package for R (Caimo and Friel, 2014) implements the methodology in this paper. It is available on the CRAN package repository at http://cran.r-project.org/web/packages/Bergm.

2 Overview of Bayesian model selection

Consider the countable model set $\mathcal{M}=\{\mathcal{M}_{1},\mathcal{M}_{2},\mathcal{M}_{3},\ldots\}$. The data $y$ are assumed to have been generated by one of the models in that set. Bayesian model selection aims at calculating the posterior model probability for model $\mathcal{M}_{m}$, $\pi(\mathcal{M}_{m}\mid y)$, where it may be of interest to obtain a-posteriori a single most probable model or a subset of likely models.

We associate each model $\mathcal{M}_{m}$ with a parameter vector $\theta_{m}$. The prior beliefs for each model are expressed through a prior distribution $p(\mathcal{M}_{m})$ $(\sum_{m\in\mathcal{M}}p(\mathcal{M}_{m})=1)$ and for the parameters within each model through $p(\theta_{m}\mid\mathcal{M}_{m})$. These specifications allow Bayesian inference to proceed by examining the posterior distribution

[TABLE]

The within-model posterior appears as $\pi(\theta_{m}\mid y,\mathcal{M}_{m})\propto f(y\mid\theta_{m},\mathcal{M}_{m})p(\theta_{m}\mid\mathcal{M}_{m})$. The constant of proportionality, termed the marginal likelihood or evidence, for model $\mathcal{M}_{m}$ is expressed by

[TABLE]

assuming a prior distribution for $\theta_{m}$ that leads to a marginal likelihood which is finite. Precise estimation of the above integral is challenging as it involves a high-dimensional integration over a usually complicated and highly variable function, so in most cases the model evidence is not analytically tractable. Knowledge of the evidence is required to deduce the posterior model probability

[TABLE]

using Bayes theorem. The probability $\pi(\mathcal{M}_{m}\mid y)$ is treated as a measure of uncertainty for model $\mathcal{M}_{m}$. Comparison of two competing models in the Bayesian setting is performed through the Bayes factor,

[TABLE]

which provides evidence in favour of model $\mathcal{M}_{m}$ compared with model $\mathcal{M}_{m^{\prime}}$. The larger $BF_{m,m^{\prime}}$ is, the greater the evidence in favor of $\mathcal{M}_{m}$ compared to $\mathcal{M}_{m^{\prime}}$. A comprehensive review of Bayes factors is presented by Kass and Raftery (1995).

In this paper we are concerned with approaches based solely on within-model simulation, where the posterior distribution within model $\mathcal{M}_{m}$ is examined separately for each $m$. Recent reviews comparing popular methods based on MCMC sampling can be found in Friel and Wyse (2012) as well as in Ardia et al. (2012).

3 Exponential random graph models

Consider the set of all possible graphs on n nodes (actors), $\mathcal{Y}$. A $n\times n$ random adjacency matrix Y on n nodes and a set of edges (relationships) describes the connectivity pattern of a graph that represents the network data. A realisation of Y is denoted with $y$ and the presence or absence of an edge (directed or undirected) between the pair of nodes $(i,j)$ is coded as

[TABLE]

An edge connecting a node to itself is not permitted so $y_{ii}=0$.

Exponential random graph models represent a general class of models for specifying the probability distribution for a set of random graphs or networks based on exponential-family theory (Wasserman and Pattison, 1996). Local structures in the form of meaningful subgraphs model the global structure of the network. ERGMs model directly the network using the likelihood function

[TABLE]

where $q(y\mid\theta)$ is the un-normalised likelihood, $s:\mathcal{Y}\rightarrow\mathbb{R}^{d}$ are sufficient statistics based on the adjacency matrix and $\theta\in\Theta\subseteq\mathbb{R}^{d}$ is the vector of model parameters (Hunter and Handcock, 2006; Snijders et al., 2006). Our focus lies on ERG models that are edge-dependent, and whose likelihood is intractable.

The evaluation of $z(\theta)$ is feasible for only trivially small graphs as this sum involves $2^{\binom{n}{2}}$ terms for undirected graphs. Recent studies on the inference of ERGMs with the Bayesian approach include Koskinen (2004), Caimo and Friel (2011), Wang and Atchade (2014), Caimo and Mira (2015), Thiemichen et al. (2016) and Bouranis et al. (2017). Bayesian model selection for exponential random graph models has been explored by Caimo and Friel (2013), Friel (2013), Thiemichen et al. (2016) and Everitt et al. (2017).

A reparameterization of (4) can express the distribution of the Bernoulli variable $Y_{ij}$ under the conditional form

[TABLE]

where $\delta_{s}(y)_{ij}=s(y^{+}_{ij})-s(y^{-}_{ij})$ denotes the vector of change statistics. The vector is associated with the dyad $y_{ij}$ corresponding to a particular pair of nodes $(i,j)$ and represents the change in the vector of network sufficient statistics when $y_{ij}$ is toggled from a 0 (no edge, $y^{-}_{ij}$) to a 1 (edge, $y^{+}_{ij}$), holding the rest of the network, $y_{-ij}=y\backslash\{y_{ij}\}$, fixed. The pseudolikelihood method, developed by Besag (1977) and applied to social networks by Strauss and Ikeda (1990), defines an approximation of the full joint distribution in (4) as the product of the full conditionals for individual observations/ dyads:

[TABLE]

where $\text{y}_{-ij}$ denotes $\text{y}\backslash\{\text{y}_{ij}\}$. The condition $i<j$ holds for undirected graphs.

4 Adjusting the pseudolikelihood

Analytical or computational intractability of the likelihood function poses a major challenge to Bayesian inference, as well as to all likelihood-based inferential approaches. A natural strategy to deal with such model intractability is to substitute the full likelihood with a surrogate composite likelihood (Lindsay, 1988; Varin et al., 2011), that shares similar properties with the full likelihood. The pseudolikelihood (Besag, 1975, 1977) is a special case of the composite likelihood and can serve as a proxy to the full likelihood when the assumption of conditional independence of the variables is reasonable.

This assumption is usually unrealistic, though. The drawback of the pseudolikelihood is that it ignores strong dependencies in the data and can, therefore, lead to a biased estimation. We propose to perform adjustments on the pseudolikelihood to obtain a reasonable approximation to the likelihood.

Adjustments to composite likelihood functions have been previously suggested by Ribatet et al. (2012). Following the proposals of Stoehr and Friel (2015) and Bouranis et al. (2017) for GRFs, we initially adjust the pseudolikelihood itself by matching its first two moments with the first two moments of the likelihood through a model-specific invertible and differentiable mapping

[TABLE]

which depends on the maximum likelihood estimate, $\hat{\theta}_{MLE}$, the maximum pseudolikelihood estimate, $\hat{\theta}_{MPLE}$, and a transformation matrix $W$. The mode and curvature-adjusted pseudolikelihood is defined as the function $y\mapsto f_{\text{PL}}(y\mid g(\theta))$.

Figure 1 displays a difference in magnitude from $f(y\mid\theta)$; a magnitude adjustment of the mode and curvature-adjusted pseudolikelihood results in the fully adjusted pseudolikelihood

[TABLE]

for some constant $C>0$. The remainder of this section provides guidelines for estimating $C$ and obtaining the mapping $g$ through the estimations of $\hat{\theta}_{MPLE}$, $\hat{\theta}_{MLE}$, and $W$.

4.1 Mode adjustment

Empirical analysis by Stoehr and Friel (2015) and Bouranis et al. (2017) showed that the Bayesian estimators resulting from using the pseudolikelihood function as a plug-in for the true likelihood function are biased and their variance can be underestimated. It is, therefore, natural to consider a correction of the mode of the pseudolikelihood approximation.

Paramount to the approach is the ability to estimate the maxima of the likelihood and the pseudolikelihood,

[TABLE]

While the MPLE is fast and straightforward to obtain because of the closed form of the pseudolikelihood, care is needed when estimating the MLE for ERGMs. We considered the Monte Carlo Maximum Likelihood Estimation (MC-MLE) procedure proposed by Geyer and Thompson (1992). Alternative procedures exist, see Hunter and Handcock (2006).

4.2 Curvature adjustment

Composite likelihoods have been previously shown to modify the correlation between the variables (Stoehr and Friel, 2015). The mapping in (5) ensures that the adjusted pseudolikelihood and the full likelihood have the same mode and aims to recover the overall geometry of the distribution (Figure 1). We choose the transformation matrix $W$ that satisfies

[TABLE]

so that the gradient and the Hessian of the log-likelihood and $\log\tilde{f}(y\mid\theta)$ are the same. It is possible to estimate the gradient and the Hessian of the log-likelihood using the following two identities:

[TABLE]

and

[TABLE]

where $\mathbb{V}_{y\mid\theta}\left[s(y)\right]$ denotes the covariance matrix of $s(y)$ with respect to $f(y\mid\theta)$. The presence of the normalising term renders exact evaluation of $\mathbb{E}_{y\mid\theta}\left[s(y)\right]$ and $\mathbb{V}_{y\mid\theta}\left[s(y)\right]$ intractable. We resort to Monte Carlo sampling from $f(y\mid\theta)$ in order to estimate these.

The Hessian matrices at the maximum of their respective function are negative definite matrices and therefore admit a Cholesky decomposition,

[TABLE]

where $M$ and $N$ are upper triangular matrices with strictly positive diagonal entries. By straightforward algebra, combining (8) and (4.2) yields $W=M^{-1}N$.

4.3 Magnitude adjustment

The magnitude adjustment aims to scale the mode and curvature-adjusted pseudolikelihood to the appropriate magnitude by performing a linear transformation of the vertical axis. The constant $C$ in (6) is defined so that $\tilde{f}(y\mid\hat{\theta}_{MLE})=f(y\mid\hat{\theta}_{MLE})$, which implies

[TABLE]

Since $z(\hat{\theta}_{MLE})$ is intractable, we unbiasedly estimate $C$ by replacing the normalising constant with an estimator which we now describe, following Friel (2013). We introduce an auxiliary variable $t\in[0,1]$ discretised as $0=t_{0}<t_{1}<\ldots<t_{L}=1$ and consider the distributions

[TABLE]

An estimate of $z(\hat{\theta}_{MLE})$ can be obtained using

[TABLE]

where $z(0)=2^{\binom{n}{2}}$ for undirected graphs and $n$ is the number of nodes. Note that in the case of a Potts/autologistic model, $z(0)=2^{N}$, where $N$ is the size of the lattice.

Importance sampling is used to estimate the ratios of normalising constants in (11). We take the un-normalised likelihood $q(y\mid t_{j}\theta)$ as an importance distribution for the "target" distribution $f(y\mid t_{j}\theta)$, noting that

[TABLE]

An unbiased importance sampling estimate of this expectation can be obtained by simulating multiple draws $y^{(j)}_{1},\dots,y^{(j)}_{K}\sim f(y\mid t_{j}\hat{\theta}_{MLE})$, yielding

[TABLE]

Increasing the number of temperatures $L$ and the number of simulated graphs will lead to a more precise estimate of $z(\hat{\theta}_{MLE})$ and will necessarily increase the computational burden. However, as we shortly illustrate, however, this does not add significantly to the overall computational cost of the adjustment procedure. We note that estimation of $z(\hat{\theta}_{MLE})$ is performed once upfront for each competing model. The estimate of $C$ is

[TABLE]

where $\hat{z}(\hat{\theta}_{MLE})$ follows from (11).

We note that estimators similar to (11) could be obtained using annealed importance sampling (Neal, 2001) or a sequential Monte Carlo algorithm (Del Moral et al., 2006). Indeed, these approaches do not require independent simulations from the likelihood for different temperatures and may therefore provide greater accuracy per computational cost and could be considered as alternatives. However, depending on the implementation of those methods, the resulting estimator may be biased.

5 Approximation of the Bayes factor

Replacing the likelihood with the unadjusted pseudolikelihood approximation in Bayes formula yields the within-model pseudo-posterior distribution

[TABLE]

In analogy to (3), the Bayes factor based on the unadjusted pseudolikelihood approximation is

[TABLE]

A naive implementation of the pseudolikelihood or any other higher-order composite likelihood is likely to give misleading marginal likelihood estimates, as we illustrate in Section 6. Having completed the adjustment steps, we propose to approximate the within-model posterior distribution by

[TABLE]

Working with (13) we can now approximate (3) by

[TABLE]

The aforementioned framework offers one possibility in the Bayesian setting to obtain an approximation to the within-model posterior distribution. Another possibility has been explored by Bouranis et al. (2017), whose ERGM experiments showed that estimation with the pseudo-posterior distribution is biased. The authors presented an algorithm to draw an approximate sample from the intractable posterior distribution $\pi(\theta\mid y)$.

Their suggested approach first samples from the pseudo-posterior distribution (12). Then an invertible and differentiable mapping $\phi:\Theta\rightarrow\Theta$ is considered to transform the entire sample $\{\theta_{i}\}^{T}_{i=1}$ so that it is a sample from an approximation of the posterior distribution, whose density is

[TABLE]

Following this approach and applying a change of variables, the model evidence $\pi(y\mid\mathcal{M}_{m})$ is approximated by

[TABLE]

and it follows that there is no gain from the transformation of the pseudo-posterior distribution when the aim is to obtain a reasonable approximation of the marginal likelihood. As such, while the correction algorithm of Bouranis et al. (2017) is appropriate for conducting Bayesian inference on the model parameters, it is not suitable for model selection.

6 Potts simulation study

The Ising model has been a popular approach to modeling spatial binary data $y=\{y_{1},\ldots,y_{N}\}\in\{-1;1\}^{N}$ on a lattice of size $N=\upsilon\times\nu$, where $\upsilon$ and $\nu$ are the number of rows (height) and columns (width) of the lattice, respectively. A lattice with $N$ nodes has $2^{N}$ possible realizations; the normalising constant $z(\theta)$ in (1) is a summation over all of the realizations and it becomes analytically unknown for moderate sized graphs. The autologistic model (Besag, 1972) extends the Ising model to allow for unequal abundances of each state value, while the Potts model (Potts, 1952) allows each lattice point to take one of $S\geq 2$ possible values/states.

In this example we investigate the efficiency of the approximation to the marginal likelihood when the likelihood is replaced by $\tilde{f}(y\mid\theta_{m},\mathcal{M}_{m})$, with a small dataset for which we can carry out exact computations. 30 realizations from an isotropic 2-state Potts model with interaction parameter $\theta=0.4$ defined on a lattice of size $15\times 15$ were exactly sampled via Friel and Rue (2007) and Stoehr et al. (2016). The sufficient statistic for the Potts model is the number of corresponding neighbors in the graph

[TABLE]

where the notation $j\sim i$ denotes that node $j$ is a neighbor of node $i$. We assume that the lattice points have been indexed from top to bottom in each column and that columns are ordered from left to right. For a first-order neighborhood model an interior point $y_{i}$ has neighbors $\{y_{i-\upsilon},y_{i-1},y_{i+1},y_{i+\upsilon}\}$; nodes situated at the boundary of the grid have less than four neighbors (Figure 2).

The model under consideration is expressed as

[TABLE]

where a diffuse Gaussian prior distribution, $\mathcal{N}\left(0,25\right)$, was assumed.

The size of the simulated lattices allows for accurate estimates of the evidence as follows: the normalising constant $z(\theta)$ can be calculated exactly with a recursive forward-backward algorithm (Reeves and Pettitt, 2004; Friel and Rue, 2007), which can then be plugged into the right hand side of (14). Numerical integration (using the trapezoidal rule) of the right hand side over a grid of $\{\theta_{i}\}_{i=1}^{M}$ values gives an accurate estimate of $\pi(y)$:

[TABLE]

This serves as a ground truth against which to compare with the corresponding estimates of the model evidence under the fully adjusted pseudolikelihood, detailed in Section 4. To reduce discretisation error we considered a sequence of length 5,000 over the interval [0, 0.8], covering the effective range of values that $\theta$ can take. The same numerical integration can be used to estimate the evidence based on the (corrected) pseudo-posterior distribution.

In the case of the Potts model the pseudolikelihood comprises the product of full-conditional distributions of each $y_{i}$,

[TABLE]

where $y_{-i}$ denotes $y\backslash\{\text{y}_{i}\}$. For a lattice of size $N$ the denominator of the pseudo-likelihood will take $2N\ll 2^{N}$ calculations to evaluate, providing a significant improvement when compared to the full likelihood. When $\theta=0$, the pseudolikelihood function is identical to the true likelihood.

In this example we estimated (11) using a ladder of 100 equally spaced path points and sampling 1500 graph statistics at each of them. The algorithm took 1.5 min to estimate the intractable normalising constant for each dataset. This setup has been empirically shown to be sufficiently accurate in this example, where the estimate of the intractable normalising constant at the mode agrees with the respective estimate using the recursive forward-backward algorithm.

In order to assess the accuracy of the approximation to the marginal likelihood based on the fully adjusted pseudolikelihood, in Figure 3 we present a scatterplot of the true value of $\log\pi(y)$ against $\log\pi_{\text{PL}}(y)$ and $\log\widetilde{\pi}(y)$ for a range of scenarios with different dependence structures. For datasets with strong dependence structure, the unadjusted pseudolikelihood is inefficient and it should be avoided for model comparisons.

On the contrary, an approximation based on the fully adjusted pseudolikelihood gives excellent performance irrespectively of the dependence structure and is virtually an exact estimate to the true evidence. We note that when we are able to calculate $z(\hat{\theta}_{MLE})$ exactly using the recursive forward-backward algorithm we have almost exact agreement between $\log\pi(y)$ and $\log\widetilde{\pi}(y)$.

7 Applications

The intractability of the ERG likelihood function leads to a doubly-intractable posterior distribution. Due to the intractable normalizing constant of the likelihood function, the estimation of the evidence is added with another layer of difficulty. Our tractable approximation to the likelihood yields an approximated posterior distribution, making it possible to use existing evidence estimation techniques in the presence of intractable multi-dimensional integrals (see Section 1).

Such within-model strategies focus on the (approximated) posterior distribution for each competing model $\mathcal{M}_{m}$ separately, with the aim to estimate their model evidence. For illustration purposes, we considered three kinds of strategies that are based on MCMC simulation: Chib and Jeliazkov’s one block Metropolis-Hastings method, thermodynamic integration (TI) and Steppingstone sampling. Below is a brief overview of these and a more detailed presentation is available at the Supplementary material. Algorithm 1 summarises the steps involved in estimating the evidence $\widetilde{\pi}(y\mid\mathcal{M}_{m})$ for each model.

• •

One block Metropolis: equation (2) is rearranged with respect to the marginal likelihood. Assuming a tractable likelihood, estimation of the marginal likelihood requires an estimate of the posterior ordinate $\pi(\theta\mid y)$, which can be found at a high-density point $\theta^{*}$ in the support of the target posterior distribution for estimation efficiency. The method is used in the case where the parameter vector $\theta$ can be updated in a single block and has been generalised to situations where the parameter vector is updated in full-conditional blocks (Chib and Jeliazkov, 2001).

• •

Power posteriors and controlled thermodynamic integration (CTI): a path sampling type method and an application of the thermodynamic integration (TI) technique from statistical physics (Friel and Pettitt, 2008). The power posterior distribution is defined to be proportional to the product of the prior and the likelihood raised to a power $t$ where $t\in[0,1]$ so that $\pi_{t}(\theta\mid y)\propto f(y\mid\theta)^{t}p(\theta)$. The inverse temperature $t\in[0,1]$ has the effect of tempering the likelihood. A sample is drawn from the power posterior distribution and it is used to estimate the intractable evidence. An improvement on the variance of the evidence estimator can be achieved through the use of control variates (Oates et al., 2016).

• •

Stepping stones sampler: uses the idea of powered posteriors, treating them as a series of intermediate distributions between the prior and the posterior. The corresponding normalising constant of the power posterior is $z(y\mid t)=\int_{\theta}f(y\mid\theta)^{t}p(\theta)\;\mathrm{d}\theta$. An estimate of the evidence, $z(y\mid t_{m}=1)$ is given by the product of $k=0,\ldots,m-1$ ratios of consecutive normalising constants, $r_{k}=z(y\mid t_{k+1})/z(y\mid t_{k})$, see Xie et al. (2011).

Computation in this paper was carried out with the statistical environment R (R Core Team, 2017) on a laptop computer with an Intel ®CoreTM i7-4500U CPU (1.80GHz) and 16GB RAM. Throughout the analysis of the network data we assumed a diffuse Multivariate Gaussian prior distribution for the model parameters, $\mathcal{MVN}\left(0_{d},100\times I_{d}\right)$, where $0_{d}$ is the null vector and $I_{d}$ is the identity matrix of size equal to the number of model dimensions, $d$, unless stated otherwise.

7.1 Benchmark algorithms

Comparisons are provided against the auto-reversible jump (Auto-RJ) exchange (Caimo and Friel, 2013) and the population exchange (Friel, 2013).

• •

The Auto-RJ is a trans-dimensional MCMC algorithm. The algorithm consists of two steps: the first (offline) step is used to sample from the posterior distribution of each competing model using the exchange algorithm (Caimo and Friel, 2011) and then to approximate the estimated posterior by Gaussian distributions determined by the first moments of each sample. The second (online) step of the algorithm makes use of the Gaussian posterior proposal estimated in the offline step as within-model proposals for the reversible jump-MCMC computation. Note that the Auto-RJ exchange requires draws from the likelihood. Since it is not possible to draw exactly from it, an auxiliary Markov chain from the tie-no-tie (TNT) sampler (Hunter et al., 2008) is used to return a draw that is approximately distributed under the true likelihood, in place of exact simulation.

• •

In Friel (2013), the Author proposes an efficient way to estimate the evidence based on the identity

[TABLE]

which holds for all $\theta^{*}\in\Theta$. Of course ${z}(\theta^{*})$ and ${\pi}(\theta^{*}\mid y)$ are unknown. Friel (2013) devises the "population" version of the exchange algorithm that yields draws, $\{\theta^{(i)}\}$, from the posterior distribution by transitioning from the prior (see Supplementary material). Additional auxiliary draws at each iteration of the MCMC scheme with the TNT sampler are used to give an estimate of the intractable normalising constant, $\hat{z}(\theta^{(i)})$, at each of these points. A kernel density approximation of the posterior, $\widehat{\pi}(\theta\mid y)$, is also found using these points. The estimator of $\pi(y)$ is obtained by averaging different estimated values of $\pi(y)$ replacing the unknown quantities in (15) by the estimated ones $\hat{z}(\theta^{(i)})$ and $\widehat{\pi}(\theta\mid y)$, for a number of draws of the posterior that are close to posterior mean.

Both applications to social networks involve comparisons against the Auto-RJ and the population exchange algorithms. The main weakness of the population exchange method is that it relies on (i) extensive simulations from the likelihood, increasing the computational burden and (ii) kernel density estimation of the target posterior distribution, which makes it impractical for use when the parameter space is high-dimensional. The application to the adolescent friendship network serves as an example where model comparisons based on unadjusted pseudolikelihoods are misleading, as opposed to comparisons based on fully adjusted pseudolikelihoods.

7.2 Gaussian random walk updates

The low-dimensional parameter space of the models considered in this paper allow to sample the target posterior distribution in one block. With the exemption of the auto-reversible jump exchange algorithm, all the other methods considered in this paper employed a random walk Metropolis strategy with a multivariate Gaussian proposal distribution.

A proposal distribution with a variance-covariance matrix in the form $\Sigma_{\lambda}=\lambda^{2}I_{d}(B_{0}+C^{-1})^{-1}$ was assumed (Chib and Jeliazkov, 2001; Martin et al., 2011), to account for possible correlations between the model parameters. Here, $\lambda\in\mathbb{R}^{+}$ denotes the Metropolis tuning scalar parameter and $B_{0}$ is the prior precision. There are two options regarding the choice of the precision matrix $C^{-1}$: (i) it is the same as the negative Hessian $-\mathcal{H}_{\text{PL}}(\hat{\theta}_{MPLE})=-\nabla^{2}_{\theta}\log{f_{\text{PL}}(y\mid\theta)}|_{\hat{\theta}_{MPLE}}$ when inference is based on the unadjusted pseudolikelihood or (ii) it is the same as the negative Hessian $-\nabla^{2}_{\theta}\log{f(y\mid\theta)}|_{\hat{\theta}_{MLE}}$ when inference is based on the adjusted pseudolikelihood.

The adjusted pseudo-posterior MCMC algorithm and the population exchange algorithm were tuned using the sample covariance matrix of the likelihood at the MLE from option (ii). Here $B_{0}$ plays the role of regularising the Hessian. This strategy assumes that the Hessian is invertible. Obviously, a bad choice of the covariance matrix $\Sigma_{\lambda}$ will have an effect on posterior inference; a short MCMC run helped us decide on a value of the Metropolis tuning parameter in order to reach a reasonable mixing rate of around 25% (Rosenthal and Roberts, 2001).

As regards the sampling algorithms based on tempered likelihoods, the parameters at all temperatures were updated jointly with Gaussian random walk proposals. It would be desirable to scale those proposals within different temperatures; it is appropriate to have wider proposals at lower temperatures so that the algorithm can explore the posterior support more effectively. When updating the parameter vector $\theta_{j}$ in temperature $t_{j}$, Friel and Wyse (2012) chose a proposal from a Gaussian distribution which was centered at $\theta_{j}$ and with standard deviation $(t^{\alpha}\tau_{p})^{-1/2}$, where $\tau_{p}$ is the proposal precision at temperature $t=1$. For their logistic regression example they chose a value of $\alpha$ such that the variability of the proposal near zero temperature would equal that of the prior. In the case of power posteriors, this was $\alpha=\log(\tau/\tau_{p})/\log(t_{1})$. This scaling led to reasonable acceptance rates in our experiments.

7.3 Karate club network

The Zachary’s Karate Club network data, displayed in Figure 4, represents a social network of friendships between 34 members of a karate club at a US university in the 1970. Three competing models were proposed by Caimo and Friel (2014) to fit the data in the presence of degeneracy:

[TABLE]

where $s_{1}(y)=\sum_{i<j}y_{ij}$ is the number of edges. The other model terms are defined below.

Shared Partnership: Let $EP_{k}(y)$, called the edgewise shared partnership statistic, denote the number of connected pairs with exactly k common neighbours. $EP_{k}(y)$ is a function of the triangle counts and as such, it is equivalent to modeling the high-order transitivities. The distribution of edgewise shared partnership can be modeled as a function of a single parameter by placing decreasing weights on the higher transitivities, leading to the geometrically weighted edgewise shared partnership (GWESP) statistic. GWESP is defined by:

[TABLE]

Geometrically Weighted Degree: Let the degree count, $D_{k}(y)$, denote the number of pairs that have exactly k common neighbours. The number of stars is a function of the degrees, therefore $D_{k}(y)$ is equivalent to modeling the k-star statistic. The geometrically weighted degree (GWD) statistic enables to model all degree distributions as a function of single parameter by placing decreasing weights on the higher degrees. GWD is defined by:

[TABLE]

The scale parameters $(\phi_{v},\phi_{u})$ specify the decreasing rates of weights placed on the higher order terms, are treated as constants and are set to $(\phi_{v}\,,\phi_{u})=(0.2\,,0.8)$. The main focus of this example will lie on the comparison between $\mathcal{M}_{1}$ and Model $\mathcal{M}_{3}$. Table 2 makes clear that there is positive evidence in favor of $\mathcal{M}_{1}$ over Model $\mathcal{M}_{3}$, implying that the effect captured by the geometrically weighted degree network statistic does not enhance the observed network.

Thirty independent MCMC experiments with the Auto-RJ exchange algorithm were run, where each simulation consisted of 500,000 iterations, discarding the initial 50,000 as part of the burn-in. The auxiliary chain consisted of 300,000 iterations. Following arguments presented in Everitt (2012), the approximate exchange algorithm (Caimo and Friel, 2011) that the Auto-RJ is using converges to the target distribution as the number of auxiliary draws tends to infinity. There are 561 possible ties in the Karate Club network and so 300,000 auxiliary iterations should ensure that the long auxiliary MCMC run from the TNT sampler returns a draw that is approximately distributed under the true likelihood.

The Auto-RJ experiments yielded an average value of $BF_{13}=13.177$ and between-model acceptance rate $\tilde{B}=5.5\%$. The low acceptance rate for the jump proposals suggests that the proposal distributions are not a good fit to each posterior model. As a result, the chain is mixing poorly and each model is not visited with the correct frequency. This, in turn, has an effect on the Bayes factor estimate.

The population exchange algorithm was implemented with 50 chains for 10,000 overall iterations after a burn-in period of 1,000 iterations, where 5,000 auxiliary iterations were used to generate an approximate draw of a network from the likelihood. A further 500 draws were used for the importance sampling estimate of the ratio of normalising constants between successive temperatures. We considered a temperature schedule $t_{i}=(i/50)^{5},\,i=0,\ldots 50$. The closest 500 MCMC draws to the posterior mean of $\theta$ were used to estimate $\pi(y)$.

Chib & Jeliazkov’s method drew a sample of 1,000,000 from the (corrected) pseudo-posterior distribution for estimating the evidence. The corresponding MCMC sampler was tuned appropriately, following the guidelines at Section 7.2 to obtain an overall acceptance rate of 20-25%.

For the tempering schemes we used a temperature ladder with 101 rungs, $t_{i}=(i/100)^{5},\,i=0,\ldots 100$. Within each temperature $t_{i}$, 30,000 samples were collected from the corresponding stationary distribution, after a burn-in period of 5,000 iterations. This extended run increased considerably the computational expense but aided towards the reduction of estimation bias (see Friel et al. (2014) for general recommendations on the number of rungs and the length of the MCMC run at each temperature $t_{i}$). We note that the additional computation related to the control variates is a negligible fraction of the total computational cost.

The importance sampling algorithm for estimating the intractable normalising constant was carried out using 500 path points. For each point the TNT sampler was run for 5,000 iterations, followed by an extra 75,000 iterations thinned by a factor of 50, yielding 1500 networks. The adjustment algorithm took 3 min for each of the two models (Table 1). The overall time spent on the adjustments of the pseudolikelihood is accounted for in Table 2.

Results, shown in Table 2, demonstrate that the magnitude and direction of the Bayes factor estimates with the adjusted pseudolikelihood agree with the Auto-RJ results. In fact, working with the unadjusted pseudolikelhood in this example gives Bayes factor estimates of similar order of magnitude, but it is not recommended to work with such an approximation as we cannot be sure about its quality beforehand for a given dataset.

Additionally, there is good agreement with the evidence estimates from the population exchange MCMC run. Our procedure gets very accurate estimates of the true Bayes factor, but in a fraction of the time, which renders it a more appealing option.

Chib and Jeliazkov’s one-block approach comes at a lower computational cost compared to TI. The Power posterior-related estimates are based on the improved trapezium rule (S.4). The potential gains made when estimating the evidence using power posteriors by correcting the numerical integration error have been illustrated previously (Friel et al., 2014). As regards the controlled thermodynamic integral estimate, the results are additionally based on a zero variance (ZV) control variates polynomial $P(\theta)$ of dimension 2 (see Supplementary material).

There is very small discrepancy between the Stepping stone estimators and the Power posterior estimators based on second-order quadrature under the long temperature ladder. CTI (degree 2) achieves a massive variance reduction in the estimator variance. We see that this variance reduction transfers to estimates of the Bayes factor themselves, where the standard deviation of the CTI estimators (degree 2) is approximately $10\times$ lower compared to estimators based on TI. Figure 5 displays that the CTI estimators (degree 2) are close to the Chib and Jeliazkov’s estimators.

7.4 Teenage Friends and Lifestyle Study

The adolescent friendship network is a subset of a friendship network collected in the "Teenage Friends and Lifestyle Study" (Pearson and Michell, 2000). The study records a network of friendships and substance use for a cohort of students in a secondary school in Glasgow, Scotland. Here we used an excerpt of 50 adolescent girls; the resulting network, displayed in Figure 6, consists of 50 nodes and 39 edges. There are four covariates associated with each node: Drug usage, Smoking status, Alcohol usage and Sport activity.

Following Caimo and Friel (2014), we restrict our attention to three node covariates: Drug usage, Sport activity, and Smoking status. We focus on the transitivity effect (expressed by the GWESP statistics with $\phi_{v}=\log{2}$), the degree heterogeneity (expressed by the GWD statistics with $\phi_{u}=0.8$) and the relationship between drug consumption and smoking (denoted by $s_{4}(y,x)=\textit{nodematch(c("smoke","drugs"))}$).

The homophily effect counts the number of edges for which two nodes share the same covariate value. When multiple relationships are studied, the "nodematch" statistic (Morris et al., 2008) counts only those on which all the covariate values match. We compare two models:

[TABLE]

where $s_{1}(y)=\sum_{i<j}y_{ij}$ is the number of edges.

Since the network is sparse, Caimo and Friel (2014) incorporated this prior knowledge by setting the parameter value for the edges statistic equal to $-1$ and set the Multivariate Gaussian prior covariance matrix of each model to be a diagonal matrix with every entry equal to 5.

The Auto-RJ was again used as a reference. The simulation was run for 400,000 iterations (equating to 50 hours of CPU time), discarding the initial 50,000 as part of the burn-in, and the auxiliary chain consisted of 200,000 iterations. The average Bayes factor (and standard deviation) based on thirty independent MCMC experiments was $BF_{21}=1.186~{}(0.015)$.

This value shows slightly positive evidence in favor of $\mathcal{M}_{2}$ (but is rather close to one), revealing that the transitivity effect and the geometrically weighted degree distribution can explain the complexity of the observed network data, being important features of the graph and that the observed network is enhanced by the homophily effect of drug usage and smoking. Similarly to Section 7.3, the Auto-RJ suffers from low acceptance rates for the jump proposals: the between-model acceptance rate was equal to $\tilde{B}=8.2\%$.

For the population exchange algorithm we considered a temperature schedule $t_{i}=(i/100)^{5},\,i=0,\ldots 100$, with 10,000 overall iterations after a burn-in period of 1,000 iterations. Each simulation consisted of an auxiliary MCMC run of length 5,000 as a proxy for an exact sampler from the likelihood and a further 500 draws were used for the importance sampling estimate of the ratio of normalising constants between successive temperatures. The closest 500 MCMC draws to the posterior mean of $\theta$ were used to estimate $\pi(y)$. This setting gave results which agree with our method, but came at a high computational cost (Table 4).

In terms of estimating the intractable normalising constant for each model based on importance sampling in our method, we used the same setting as in Section 7.3. In this example the adjustment algorithm took 5.5 min for each of the two models (Table 3). The overall time spent on the adjustments of the pseudolikelihood is accounted for in Table 4.

With the Chib & Jeliazkov’s method 1,000,000 MCMC updates from the (corrected) pseudo-posterior distribution were run. Regarding the TI-related schemes, we ran 30,000 Metropolis-Hastings updates for each temperature, after a burn-in period of 5,000 iterations, while the tempering scale was partitioned as $t_{i}=(i/100)^{5},\,i=0,\ldots 100$. In total, 3,030,000 iterations were used to estimate the evidence for each model.

For this example we observe that the Bayes factor estimates based on the adjusted pseudolikelihood agree with the Auto-RJ results (Table 4). There is also a good agreement between the Stepping stone estimators and the Power posterior estimators before applying the ZV scheme. In general, all estimators exhibit internal consistency; a 20-fold reduction in the evidence estimator variance is observed when applying the ZV control variates scheme.

The Teenage Friends and Lifestyle Study network offers a scenario where any model comparisons based on the unadjusted pseudolikelihood are quite misleading, as the Bayes factor based on the unadjusted pseudolikelihood is $<1$, while the Bayes factor based on the adjusted pseudolikelihood is $>1$. There is a striking difference between the unadjusted pseudolikelihood-based estimates of the evidence and those estimates based on the fully adjusted pseudolikelihood. All Bayes factor estimates slightly favour $\mathcal{M}_{1}$, which is in contrast with the Auto-RJ results and the results based on fully adjusted pseudolikelihoods. Therefore, we can conclude that the benefits of conducting model selection for this network based on the unadjusted pseudolikelihood approximation are reduced. All in all, we highly recommend approximation of the evidence with the fully adjusted pseudolikelihood, which comes at a negligible computational cost.

8 Discussion

In this paper we have presented a novel approach to marginal likelihood estimation of models with intractable normalising constants, which we applied to the challenging setting of exponential random graph models for social network analysis. We approximated a doubly-intractable posterior distribution with an intractable likelihood by a "standard" singly intractable posterior distribution with a tractable likelihood approximation. Our methodology is highly compatible with a plethora of evidence estimation techniques from the Bayesian toolbox.

Our experiments suggest that the one-block Metropolis-Hastings approach yields marginal likelihood estimators with low variability. It comes at a low computational cost and is suitable because it requires no further tuning of the MCMC algorithm and the low-dimensional parameter spaces in our experiments allow the parameter vector to be updated in a single block. For higher-dimensional MCMC problems, though, multi-block Metropolis-within-Gibbs updating strategies will be more suitable (Chib and Jeliazkov, 2001). The Power posteriors algorithm can come at a higher computational cost, depending on the inverse temperature scheme. We suggest the use of the improved trapezoidal scheme of Friel et al. (2014) and the control variate technique (Oates et al., 2016) to vastly improve the statistical efficiency of the evidence estimate. The Stepping stones estimators suffer from high variability.

We note that our approach avoids the heavy computational burden of repeated likelihood simulations, as in Friel (2013) and Everitt et al. (2017). Here the likelihood simulations are performed only once for each competing model. The empirical results presented above suggest that our approach gives similar estimates of the Bayes factor to the computationally intensive population exchange algorithm, but at a fraction of the time. Simulation procedures for approximating a solution to the likelihood equation are more challenging when larger datasets are considered, but any simulation approach with increased dependence on likelihood simulations will be infeasible under these conditions. Our work should offer a more scalable approach.

Overall, the reasonable computational effort that needs to be put when working with our likelihood approximation extends the applicability of the proposed approach to other complex models like Gaussian Markov random fields and autologistic models; it will be particularly interesting to implement our proposed framework to large grids and networks with hundreds of nodes or models with many parameters.

SUPPLEMENTARY MATERIAL

Evidence estimation strategies:

Detailed presentation of the strategies used in this paper for estimating the marginal likelihood, based on MCMC simulation. (.pdf file)

R code:

R programs that can be used to replicate the Potts model study and exponential random graph example in Section 7 of this article. Please see the file README.txt contained within the accompanying zip file for more details.