Introduction

This is a short software tutorial accompanying the paper “A Stochastic Block Prior for Clustering in Graphical Models”. The paper introduces the Stochastic Block Model as a prior on the network structure of the ordinal Markov random field graphical model (Marsman et al., 2025). In addition to offering an approach for incorporating the theoretical assumption of clustering directly into the statistical model, the method provides a principled way to infer the number of clusters in the network given the data, and to test hypotheses about clustering. Although the method presented in the paper is limited to the ordinal Markov random field, the software implementation has since been extended to Gaussian graphical models and graphical models for mixtures of continuous and ordinal data.

In this brief tutorial, an example analysis is first illustrated using the R package bgms (Marsman et al., 2023). The same analysis is then demonstrated using the user-friendly wrapper package easybgm (Huth et al., 2024). The package easybgm provides a more accessible interface for conducting the analysis, particularly for users without a strong background in R. These functionalities are also available in the newest Network module of the open-source statistical software JASP. Please install the Beta version of the module. Then, under Priors, select the Stochastic Block Model, and additional clustering output options will appear.

Suppose the aim is to evaluate whether the network structure estimated from empirical data exhibits clustering. Within the Bayesian graphical modeling framework, this can be examined by estimating a representative node allocation vector from the posterior distribution. This allocation summarizes the most plausible cluster membership of each node given the data and thereby provides a direct indication of potential community structure in the network. In addition, the model yields posterior probabilities over the number of available components, denoted by \(B\). The Bayes factors below test hypotheses about \(B\). The node allocations and co-clustering matrix describe how the observed nodes are grouped.

The first step is to estimate the model by passing a stochastic block model (SBM) prior object to the edge_prior argument, using the sbm_prior() function. When using the SBM as a prior on the network structure, several of its arguments are important:

• Shape hyperparameters for the Beta prior distributions set on the probabilities for the inclusion of within and between block edges
  • alpha and beta (within clusters)
  • alpha_between and beta_between (between clusters)
• These hyperparameters govern the expected density of edges within and between clusters.
• Larger values of alpha and smaller values of beta encourage denser connections within clusters.
• Conversely, smaller values of alpha_between and larger values of beta_between encourage sparser connections between clusters.
• To calculate the proportion of connected edges implied by a specific set of hyperparameter values (for the within or between edges), we can use the expression for the expected value of the beta distribution:

\[\frac{\alpha}{\alpha + \beta}.\]

• If the substantive hypothesis is that connections are denser within clusters than between clusters, choose hyperparameters that express this expectation and examine sensitivity to those choices. Setting all four hyperparameters to 1 gives uniform priors on both inclusion probabilities and does not express a preference for denser within-cluster connections.

• The rate hyperparameter of the shifted Poisson prior on the number of available components

  • lambda
    

  • The prior is \(B - 1 \sim \text{Poisson}(\lambda)\), so the expected number of available components is exactly \(1 + \lambda\), and the prior probability of one available component is \(p(B = 1) = e^{-\lambda}\). For example:

    lambda	Expected available components, \(E[B]\)	\(p(B = 1)\)
    0.5	1.5	0.61
    1	2	0.37
    2	3	0.14
    3	4	0.05
    4	5	0.02

  • For the positive values of \(\lambda\) accepted by the software, an expected number of exactly one available component is only approached as \(\lambda \to 0\). A small value such as lambda = 0.5 could be reasonable when one wishes to favour \(B = 1\) a priori while still allowing more components: it assigns approximately 61% prior probability to \(B = 1\). This is a statement about prior preference, not a guarantee of more conservative Bayes-factor decisions, since Bayes factors also adjust for prior odds.

An adjusted count prior. The paper describes a zero-truncated Poisson prior, whereas the new bgms 0.2.0.0 uses a shifted Poisson prior. Both allow one or more available components, but the same lambda implies different expected counts. For example, lambda = 1 gives an expected count of approximately 1.58 under the paper’s prior and 2 under the shifted prior. With the shifted prior, the expected count is simply \(1 + \lambda\): use lambda = 1, 2, or 3 for expectations of 2, 3, or 4 available components, respectively. Use the table above when choosing lambda for new analyses and sensitivity checks when using bgms version 0.2.0.0 or higher.

• The concentration parameter of the Dirichlet distribution on the probabilities of allocating the nodes into clusters
  • dirichlet_alpha
  • Conditional on \(B\), this hyperparameter controls the distribution of allocation weights and thereby influences occupied cluster sizes.
    • Larger values concentrate the weights around equal shares, encouraging more balanced allocations.
      
    • Smaller values favour uneven weights, making empty components and unequal occupied cluster sizes more likely.
      
    • The default value of 1 gives a uniform distribution over the allocation-weight simplex. It does not make all partitions of the nodes equally likely.

It is strongly recommended that researchers carefully examine their choice of hyperparameter values when conducting analyses and perform prior sensitivity checks. For further discussion and practical guidance, see, for example, Marsman et al. (2023), Huth et al. (2024) and Sekulovski et al. (2024)

In the code examples below, data is a placeholder for the actual dataset to be analyzed. This dataset should be an \(n \times p\) matrix or data frame, where \(n\) denotes the number of observations and \(p\) the number of variables to be included in the network.

Suppose we have a dataset with 200 observations and 10 variables, and theoretical considerations suggest a small number of clusters with denser connections within clusters than between them. An illustrative prior specification is:

• alpha = 8,
• beta = 1,
• alpha_between = 1,
• beta_between = 8,
• lambda = 1,
• dirichlet_alpha = 1

This prior has an expected 2 available components. With 10 nodes and dirichlet_alpha = 1, the expected occupied count is approximately 1.755, rather than exactly 2. It also expresses an expectation of dense within-cluster connections (8 / (8 + 1) = 0.89) and sparse between-cluster connections (1 / (1 + 8) = 0.11).

Clustering Analysis using the package bgms

The bgms examples and count-summary discussion below refer to version 0.2.0.0. Install bgms and record the version used for the analysis.

install.packages("bgms")
library(bgms)
packageVersion("bgms")
?bgm # for more details on the function and its arguments

# run the model
fit <- bgm(data,
           edge_prior = sbm_prior(alpha           = 8,
                                  beta            = 1,
                                  alpha_between   = 1,
                                  beta_between    = 8,
                                  lambda          = 1,
                                  dirichlet_alpha = 1),
           update_method = "adaptive-metropolis")

Before interpreting any of the output, it is important to check that the sampler has converged. Please see the package help files for more details.

Now, first we check the estimated posterior cluster allocations of the nodes. The output is a vector whose length equals the number of variables (columns) in the dataset. There are two possible options: the posterior mean vector and the posterior mode vector. To examine this, we can inspect the allocation vectors, which we obtain with the extract_sbm() function.

sbm <- extract_sbm(fit)

cat("Posterior mean:\n")
print(sbm$posterior_mean_allocations)

cat("\nPosterior mode:\n")
print(sbm$posterior_mode_allocations)

We can also inspect the reported probabilities for the available-component count \(B\). In bgms 0.2.0.0, this summary has rows for \(B = 1, \ldots, p\). It is not the distribution of the occupied count \(T\): for example, \(B = p\) does not imply that each variable forms its own cluster, because some components can be empty.

sbm$posterior_num_blocks

The Bayes factors below use the corrected easybgm::clusterBayesfactor() helper, which also accepts a raw bgms fit and reconstructs the count probabilities under the full shifted Poisson prior to numerical accuracy. This implementation is available in the GitHub version. It is designed for the prior and summary convention of bgms 0.2.0.0; do not apply it indiscriminately to historical fits. Install an easybgm version containing this correction and record the package versions with the analysis:

install.packages("remotes")
remotes::install_github("KarolineHuth/easybgm")
packageVersion("easybgm")
sessionInfo()

The clustering Bayes factors are defined as follows:

\[\text{BF}_{10} = \underbrace{\frac{p(\mathcal{H}_1\mid \text{data})}{p(\mathcal{H}_0 \mid \text{data})}}_{\text{Posterior Odds}} \bigg/ \underbrace{\frac{p(\mathcal{H}_1)}{p(\mathcal{H}_0)}}_{\text{Prior Odds}}\]

Depending on the hypotheses, we have two types of Bayes factors:

• The Bayes factor comparing \(\mathcal{H}_1: B > 1\) with \(\mathcal{H}_0: B = 1\). The null implies that all observed nodes occupy one cluster. The alternative allows multiple occupied clusters, but also allows one occupied cluster with additional empty components.

The posterior odds under the full count prior are

\[\frac{p(B > 1 \mid \text{data})}{p(B = 1 \mid \text{data})} = \frac{1-p(B = 1 \mid \text{data})}{p(B = 1 \mid \text{data})}.\]

After calculating the posterior odds, we need to divide by the prior odds in order to obtain the Bayes factor. Under the shifted Poisson prior we have \(p(B = 1) = e^{-\lambda}\), so the prior odds for this Bayes factor are defined as:

\[\frac{1 - e^{-\lambda}}{e^{-\lambda}} = \exp(\lambda) - 1.\]

Dividing the posterior odds by these prior odds gives \(\text{BF}_{10}\). The corrected helper obtains the fitted hyperparameters and computes both odds internally, including the correction for the full count support:

BF_10 <- easybgm::clusterBayesfactor(fit, type = "complement")
cat("Bayes factor in favor of H1:\n")
BF_10  
cat("Bayes factor in favor of H0:\n")
1/BF_10 # inverse of the Bayes factor in favor of H1

• The Bayes factor comparing two specific cluster counts, \(\mathcal{H}_1: B = b_1\) and \(\mathcal{H}_2: B = b_2\). For example, \(b_1 = 1\) and \(b_2 = 2\). The posterior odds are \(p(B=b_1 \mid \text{data})/p(B=b_2 \mid \text{data})\).

The prior odds for this Bayes factor are given by:

\[\lambda^{b_1 - b_2} \cdot \frac{(b_2 - 1)!}{(b_1 - 1)!}\]

Again, the Bayes factor is the posterior odds divided by the prior odds. Use the corrected helper to calculate it:

b1 <- 1
b2 <- 2
BF_12 <- easybgm::clusterBayesfactor(fit, type = "point", b1 = b1, b2 = b2)
cat("Bayes factor in favor of H1:\n")
BF_12  
cat("Bayes factor in favor of H2:\n")
1/BF_12 

Researchers can additionally examine the posterior co-clustering matrix, a symmetric matrix that contains, for each pair of variables, the proportion of posterior samples in which the two variables are assigned to the same cluster. This matrix offers an informative visual summary of clustering uncertainty. Diffuse or heterogeneous patterns suggest that certain variables frequently change cluster membership or do not have a well-defined allocation.

install.packages("pheatmap")
library(pheatmap)
pheatmap(sbm$posterior_mean_coclustering_matrix, 
         cluster_rows = FALSE, 
         cluster_cols = FALSE,
         color = colorRampPalette(c("#B8860B", "#009E73"))(100))

Clustering Analysis using the package easybgm

The fitting and summary steps can also be performed through easybgm. Use the GitHub version with the corrected count Bayes-factor helper described above:

install.packages("remotes")
remotes::install_github("KarolineHuth/easybgm")
library(easybgm)
library(bgms)

fit <- easybgm(data, 
               type = "ordinal",
               edge_prior = sbm_prior(alpha           = 8,
                                      beta            = 1,
                                      alpha_between   = 1,
                                      beta_between    = 8,
                                      lambda          = 1,
                                      dirichlet_alpha = 1),
               update_method = "adaptive-metropolis")

The summary function displays the main results and calculates the Bayes factor automatically. Convergence and prior sensitivity still need to be assessed before interpreting it.

summary(fit)

Following the edge overview (for more details, see Huth et al. (2024)), the summary provides the reported available-component probabilities, the estimated node memberships, and the Bayes factor comparing \(B > 1\) with \(B = 1\).

In case researchers wish to calculate the second type of Bayes factor, similar to the steps above, they can simply use the function clusterBayesfactor, by setting the argument type = "point" and specifying the values for \(b_1\) and \(b_2\).

easybgm::clusterBayesfactor(fit, type = "point", b1 = 1, b2 = 2)

The posterior coclustering matrix can be plotted in a similar way as shown above for bgms.

library(pheatmap)
pheatmap(fit$sbm$posterior_mean_coclustering_matrix, 
         cluster_rows = FALSE, 
         cluster_cols = FALSE,
         color = colorRampPalette(c("#B8860B", "#009E73"))(100))