bgmCompare Internals

How bgmCompare() works under the hood: multi-group posterior, contrast parameterization, and sampler implementation. Source: src/models/bgmCompare/ directory.

Differences from single-group bgm

The bgmCompare engine shares the same OMRF likelihood (pseudolikelihood over ordinal/Blume-Capel variables) and the same NUTS/Metropolis sampler options, but differs in three key ways:

Parameter decomposition — Group-specific parameters $\boldsymbol{\omega}^{(g)}$ are decomposed into a shared baseline $\boldsymbol{\omega}$ and contrast parameters $\boldsymbol{\delta}^{(k)}$, linked through a contrast matrix.
Separate sampler — The bgmCompare engine uses its own chain runner (run_gibbs_sampler_bgmCompare()) rather than the general run_mcmc_chain(), because the multi-group structure requires group-specific likelihood evaluations and contrast selection. The gibbs in the runner’s name is historical: it hosts the same NUTS and adaptive-Metropolis updates this page describes; what is bespoke is the chain loop around them, not the samplers inside.
Contrast selection — Instead of edge selection (which edges exist), bgmCompare performs contrast selection: for each parameter, the baseline component is always present, and the spike-and-slab mechanism decides whether contrast parameters are active or collapsed to zero.

Contrast-matrix parameterization

For $G$ groups, the per-group parameter vector $\boldsymbol{\omega}^{(g)}$ is expressed as:

\[ \omega_{ij}^{(g)} = \omega_{ij} + \sum_{k=1}^{G-1} c_{gk}\, \delta_{ij}^{(k)} \]

where $\mathbf{C}$ is a $G \times (G-1)$ projection (contrast) matrix. The baseline $\omega_{ij}$ captures the shared structure; the contrast parameters $\delta_{ij}^{(k)}$ are the projections onto the $G-1$ orthogonal contrast directions. For two groups, the single contrast reduces to a scaled group difference. Category thresholds are decomposed in the same way.

The flat parameter vector (used by NUTS) groups by component type, not per effect:

[ overall_mains | overall_pairs | active_main_diffs | active_pair_diffs ]

Contrast parameters appear only when their inclusion indicator is 1; inactive contrasts are excluded from the vector entirely.

Multi-group log-posterior and gradient

The log-pseudoposterior sums over groups. For each group $g$:

Reconstruct group-specific parameters via the contrast matrix
Compute the pseudolikelihood using group $g$’s data (sufficient statistics: category counts, Blume-Capel stats, pairwise stats)
Evaluate the gradient of the pseudolikelihood with respect to the active parameters

The gradient() function combines three components:

Observed-data gradient (gradient_observed_active()) — projects group-level sufficient statistics into the active-parameter layout using the contrast matrix. Precomputed from counts_per_category, blume_capel_stats, and pairwise_stats.
Expected-data gradient — computed inline within gradient() by building group-specific parameters via the contrast projection, evaluating softmax probabilities, and subtracting expected counts.
Prior gradient — penalties on both baseline and contrast pairwise parameters through the polymorphic prior objects: the baseline family comes from interaction_prior (normal_prior(scale = 1) by default) and the contrast family from difference_family ("Normal" by default, or "Cauchy") at scale difference_scale. See Parameter Priors.

Sampler implementation

run_gibbs_sampler_bgmCompare() orchestrates one MCMC chain. Each iteration:

Missing data imputation (if na_action = "impute") — Impute missing observations per group using the current parameters. Updates the observation matrix, category counts, Blume-Capel statistics, and pairwise statistics.
Contrast selection (if active) — A Metropolis step decides whether each contrast parameter is included or collapsed to zero (update_indicator_differences_metropolis_bgmcompare()). Each decision also records a Rao-Blackwellized inclusion draw and feeds the per-indicator odds accumulators, as in the single-group samplers (see Rao-Blackwellized inclusion draws).
Parameter updates — With adaptive Metropolis, main effects and pairwise effects are updated separately (update_main_effects_metropolis_bgmcompare(), then update_pairwise_effects_metropolis_bgmcompare()). With NUTS, all active parameters are updated jointly via the vectorized log-pseudoposterior.
Difference-prior update (when difference_selection = TRUE) — Branches on the prior family:
- Beta-Bernoulli: pooled scalar Beta conjugate update on the shared inclusion probability.
- Stochastic-Block: pairwise off-diagonal updates delegated to BaseEdgePrior::update() (collapsed Gibbs over allocations and block probabilities); when main_difference_selection = TRUE, a separate Beta-Bernoulli conjugate update is run on the diagonal using the SBM’s within-cluster alpha and beta.
- Bernoulli: no-op (fixed inclusion probability).
The difference-prior updates are always the plain (uncorrected) conjugate forms: the normalizing-constant corrections of the continuous models do not apply here because there is no precision matrix, hence no determinant tilt (see SBM Prior).

The warmup schedule follows the same multi-stage structure as the single-group model (see Warmup Schedule).

Sweep state and normalizer cache

The Metropolis sweeps share a persistent CompareSweepState (bgmCompare_state.h) holding three things: the per-group observation matrices, the group-specific effective pairwise weights, and the rest-score matrices. The state is kept consistent incrementally. An accepted pairwise or indicator move updates the weights and two residual columns per group; an imputation or NUTS update triggers a full weight rebuild.

The state also caches the per-variable, per-group log-normalizer sums of the pseudolikelihood, under three rules. A variable’s row stays valid while its main effects and residual columns are unchanged. An accepted move stores the proposed-state normalizers. An indicator flip invalidates only the affected variables. Together this turns the per-move cost of a sweep from a full normalizer recomputation into two column updates and a cache read.

Active parameter tracking

Because contrast selection can turn parameters on and off, the total_length() function recomputes the active parameter dimension at each iteration. Index maps (build_index_maps()) provide the mapping between the flat active-parameter vector and the structured (variable, group) parameter layout.

When a contrast parameter is “off” (indicator = 0), it is excluded from the NUTS block. The step size and mass matrix adapt to the current active dimension.

Output handling

The bgmCompareOutput struct stores:

Field	Content
`main_samples`	Main-effect posterior samples (iterations $\times$ thresholds $\times$ groups)
`pairwise_samples`	Pairwise-effect posterior samples (iterations $\times$ pairs $\times$ groups)
`indicator_samples`	Contrast inclusion indicators (iterations $\times$ parameters)
`rb_inclusion_samples`	Rao-Blackwellized inclusion draws, same layout as `indicator_samples`; NaN marks indicators that were not selected
`rb_counts`	Per-indicator odds accumulators (parameters $\times$ 4) on the alpha scale: `[n01, n10, n0_visits, n1_visits]`
`allocation_samples`	Cluster allocations (iterations $\times$ variables; SBM difference prior only)
`treedepth_samples`	Tree depth per iteration (NUTS only)
`divergent_samples`	Divergent transition flags (NUTS only)
`energy_samples`	Hamiltonian energy per iteration (NUTS only)
`accept_prob_samples`	Mean per-trajectory acceptance probability (NUTS only)

As in the single-group ChainResult, buffers are preallocated with “not sampled” sentinels (−1 for integer fields, NaN for floating ones), so interrupted runs return recognizable placeholders.

The R output builder (build_output_compare()) splits the posterior means into:

$group_posterior_means — Posterior means of the group-level parameters $\boldsymbol{\omega}^{(g)}$, reconstructed from the baseline $\omega_{ij}$ and contrast parameters $\delta_{ij}^{(k)}$
$difference_posterior_means — Posterior means of the contrast parameters $\delta_{ij}^{(k)}$
$difference_inclusion_probabilities — Proportion of iterations where each contrast indicator was included