bgms

**Pairwise parameters are reported on the association scale.**

• bgm() stores and reports pairwise interaction parameters on the association scale: the coefficient omega that enters every conditional as 2 * omega * x. In 0.1.6.3 the same coupling was stored as sigma = 2 * omega. Reported pairwise effects — coef(), summary(), extract_pairwise_interactions(), and the raw posterior samples — are therefore about half their 0.1.6.3 values for the same data and model. Main effects, indicators, and inclusion probabilities are unaffected in scale. Code that reads raw pairwise samples, or that compares numbers against a 0.1.6.3 run, needs the factor of two applied; there is no automatic conversion.
• bgmCompare() moves to the association scale with bgm(), in the same release and by the same factor: the two entry points were mutually consistent in 0.1.6.3 and remain so. Pairwise effects reported by coef(), summary(), extract_pairwise_interactions(), extract_group_params(), and the raw posterior samples are therefore about half their 0.1.6.3 values, exactly as for bgm(). **A stored bgmCompare() fit from an earlier version should be refit rather than rescaled**, because the prior acted on the old, doubled coordinate: the same difference_scale or interaction_prior number is now a tighter prior than it was, and difference inclusion Bayes factors move accordingly. The calibration of the difference_scale default under the association-scale parameterization is still under study, so difference verdicts near a decision threshold should be read as scale-contingent.
• simulate_mrf() reads its pairwise argument on the association scale too: the discrete Gibbs conditional now accumulates 2 * (score - baseline) * pairwise where 0.1.6.3 accumulated (score - baseline) * pairwise. **A pairwise matrix that a 0.1.6.3 script passed in now means twice the coupling it used to**, so simulated data are more strongly dependent than before at the same input. Halve a hand-specified pairwise matrix to reproduce 0.1.6.3 output; a matrix taken from a 0.2.0.0 bgm() fit needs no adjustment, because the fit is on the same scale as the simulator. The continuous (variable_type = "continuous") branch takes a precision matrix and is unaffected.

**Changed defaults.**

• The default interaction prior is normal_prior(scale = 1) in bgm(), where 0.1.6.3 used a Cauchy at pairwise_scale = 2.5; sample_ggm_prior() uses the same default, so the GGM prior chain matches bgm(). bgmCompare() takes the same default for its baseline pairwise interactions, so a compare fit and separate bgm() fits at stated defaults now price a baseline interaction alike. Pass interaction_prior = cauchy_prior(scale = 2.5) to bgm() for the 0.1.6.3 prior, remembering that the scale now acts on the association coordinate. Under the joint precision-graph specification the normalizer correction table is keyed on the interaction prior, so the first fit at the new default cell (with a Beta-Bernoulli or Stochastic-Block edge prior) builds and caches a fresh table once.
• bgmCompare() now defaults to interaction_prior = normal_prior(scale = 1) for the baseline pairwise interactions and to difference_family = "Normal" for the group differences, both of which were Cauchy in 0.1.6.3 and through 0.2.0.0's development; the baseline default now matches bgm(), so the two entry points no longer ship different priors for an identically named argument. Results under the defaults change: a fit at the new defaults is a different model from a fit at the old ones, and difference verdicts on categories one group never observed move the most, because a Normal slab bounds an unidentified threshold difference far more tightly than a Cauchy one. cauchy_prior(scale = 1) and difference_family = "Cauchy" remain fully available and reproduce the previous behaviour.
• iter and warmup default to 2e3 in bgm() and bgmCompare(), up from 1e3. A default call therefore runs about twice as long and returns twice as many draws as it did in 0.1.6.3. Set iter and warmup explicitly to recover the old budget.
• The NUTS target_accept default is 0.80 in both bgm() (was 0.60) and bgmCompare() (was 0.65). Higher target acceptance means smaller steps: expect more leapfrog steps and a longer run per iteration, in exchange for fewer divergences. The adaptive-Metropolis default is unchanged at 0.44. Pass target_accept to restore the old value.
• On continuous and mixed models the prior on the precision diagonal defaults to exponential_prior(eta = 1), where eta is the rate in the frame in which the pairwise slab has unit scale; the raw rate is derived at fit time as eta / s for interaction-prior scale s. At the default unit slab scale this is the same distribution as gamma_prior(shape = 1, rate = 1), but with a non-unit slab scale the diagonal rate now co-scales. One consequence: a scale-free interaction prior (beta_prime_prior()) on a continuous model requires the diagonal prior to be given with rate rather than eta.

**Removed and reshaped arguments.**

• update_method = "hamiltonian-mc" has been removed; use update_method = "nuts". NUTS adapts trajectory length dynamically and is more reliable, especially with edge selection on GGM models. The hmc_num_leapfrogs argument has been removed with it. A call passing either now fails. On continuous data the new conjugate update_method = "gibbs" (see New features) is a second option that needs no step-size tuning at all.
• The standardize argument of bgm() and bgmCompare() no longer does anything: pairwise interactions are on the association scale and share one prior scale, so the per-pair adjustment by the product of the two variables' maximum scores has been dropped. The argument stays as a deprecated formal so the failure is informative rather than an unused-argument error. standardize = FALSE, the old default, warns and proceeds, because it is what the sampler already does; standardize = TRUE errors and points to setting the scale directly through interaction_prior (and difference_scale for bgmCompare()). There is no per-pair equivalent, and none is planned: standardization belongs in the model rather than in an argument, and a forthcoming g-prior formulation of the interaction prior addresses it there. The per-pair maximum-score adjustment is retired in favour of that direction.
• Priors are now objects. The scalar arguments pairwise_scale, main_alpha, main_beta, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, and lambda are deprecated formals: passing one warns and is translated into the corresponding prior object, so old calls keep working for now. edge_prior takes a prior object (bernoulli_prior(), beta_bernoulli_prior(), sbm_prior()); the old character values are still accepted through a deprecation branch. The replacements are interaction_prior, threshold_prior, means_prior, precision_scale_prior, and edge_prior. Note that pairwise_scale = s translates to cauchy_prior(scale = s), preserving the 0.1.6.3 family, so a deprecated call and a call at the new default are not the same model.
• extract_category_thresholds() is deprecated in favour of extract_main_effects(), which covers category thresholds, continuous means, and precision diagonal entries in one call.

**Changed output and semantics.**

• extract_ess() reports the Rao-Blackwellized effective sample size for the edge (or difference) indicators — the n_eff column of the fit summary's inclusion table — where 0.1.6.3 returned the indicator chain's transition-based ESS. A 0.1.6.3-era call therefore gets different numbers with no warning. The extractor family is now coherent: the Rao-Blackwellized inclusion probability (extract_posterior_inclusion_probabilities()), its R-hat (extract_rhat()), and its ESS all describe the same estimate. Fits made with bgms < 0.2.0.0 carry no Rao-Blackwellized draws, so a default call on such a fit falls back to the transition ESS without warning, while an explicit estimator = "rb" errors.
• The Rao-Blackwellized estimate is the canonical inclusion probability throughout a fit: extract_posterior_inclusion_probabilities() defaults to estimator = "rb" (with "raw" still available) and fit$posterior_mean_indicator reports it. Numbers quoted from a 0.1.6.3 fit will not reproduce exactly; the two estimate the same quantity, the Rao-Blackwellized one with lower variance.
• The fit summary's inclusion table has different columns. It now reports the Rao-Blackwellized inclusion probability with a continuous Monte Carlo standard error (mcse), effective sample size (n_eff), and split-R-hat (Rhat), beside the per-direction transition counts (n0->1, n1->0). The transition-based n_eff_mixt column is gone (see Deprecated). The three Rao-Blackwellized columns are NA where the Rao-Blackwellized draws are constant to double precision, which places the inclusion probability at its numerical bound, and are computed everywhere else.
• On a bgmCompare() fit the Rao-Blackwellized inclusion probabilities are NA for indicators the sampler never updated — the main-effect differences, unless main_difference_selection = TRUE. **A previously all-numeric matrix can now contain NA**, which breaks downstream code that assumes otherwise. The NA is honest, and the contrast with bgm() is the point: an always-included edge indicator in bgm() is still *proposed* at every iteration, so the Rao-Blackwellized accumulators see its conditional inclusion odds and report a finite Bayes factor. An unselected main-difference indicator is never proposed at all, so there is no Rao-Blackwellized quantity to report — the entry is NA exactly when the indicator was never visited, not because a number was lost. Use estimator = "raw" for the old all-numeric behaviour, or filter on is.na().
• R-hat is the classic split-R-hat (Gelman et al. 2013 / Stan; Rhat = sqrt(var_plus / W) with var_plus = (n-1)/n * W + B/n). 0.1.6.3 applied the Brooks-Gelman degrees-of-freedom adjustment used by coda::gelman.diag, which collapsed to a data-independent sqrt(5/3) ~ 1.29 on nearly-saturated binary edge (or difference) indicators — the typical shape of a decisive edge, where one split sub-chain carries a brief excursion and the rest are constant — so the most decisive edges were flagged as unconverged. Reported R-hat changes slightly for all parameters (the adjustment factor was near 1 for healthy chains) and substantially for near-saturated indicators, which now report values near 1, or NA when every chain is identical. Chains stuck constant at different values report +Inf instead of NA. See the diagnostics vignette for how to read indicator R-hat alongside the transition counts and the inclusion probability.
• Fitted objects from bgm() and bgmCompare() are S7 class objects (new dependency: S7). All existing $, [[, and names() access patterns continue to work. When an incompatible easybgm version is loaded, bgms returns plain S3 lists for backwards compatibility; this shim will be removed in a future release. The compatibility mode is chosen when the fit object is constructed, so load easybgm before fitting when the fit is meant for use with an older easybgm.
• summary() on a bgms fit returns an additional quadratic element (the precision diagonal, reported on the residual-variance scale for continuous and mixed fits) and carries a main_label attribute.

New features

**New model classes.**

• Gaussian graphical models: bgm(x, variable_type = "continuous") fits a GGM with Bayesian edge selection. Sampling uses NUTS on a free-element Cholesky (theta-space) parameterization of the precision matrix, which keeps the precision matrix positive-definite by construction; adaptive-metropolis is also available.
• A conjugate GGM Gibbs sampler: update_method = "gibbs" fits a GGM with a row-by-row update of the precision matrix, with or without edge selection. It needs no step-size or proposal tuning and supports a Normal or Cauchy prior on the edges. Continuous data only.
• With update_method = "gibbs" and edge selection, the first 15% of warmup runs the full model so the precision matrix settles, and edge selection is active for the remaining 85%. Both windows scale with the warmup budget.
• Mixed MRF models: bgm() accepts a per-variable variable_type vector that mixes "ordinal", "blume-capel", and "continuous" types, to estimate networks with both discrete and continuous variables. simulate.bgms() and predict.bgms() support mixed models.
• Missing-data imputation now covers the new model classes: na_action = "impute" integrates over missing values during MCMC sampling for continuous and mixed models as well as ordinal ones.

**Prior specification.**

• Priors are constructed by dedicated functions and passed as objects: normal_prior(), cauchy_prior(), and beta_prime_prior() for parameter priors; gamma_prior() and exponential_prior() for the precision diagonal; bernoulli_prior(), beta_bernoulli_prior(), and sbm_prior() for the edge prior. Each has a print() method. The corresponding bgm() formals are interaction_prior (pairwise slab), threshold_prior (category thresholds, default beta_prime_prior(0.5, 0.5) — the 0.1.6.3 main_alpha/main_beta pair), means_prior (continuous means, default normal_prior(scale = 1)), precision_scale_prior (precision diagonal), and edge_prior.
• gamma_prior() and exponential_prior() accept the rate as eta, the rate on the precision diagonal in the frame where the pairwise slab prior has unit scale. The raw rate is derived at fit time as eta / s for interaction-prior scale s, so the prior geometry is held fixed when the slab scale changes; at fixed eta, graph and partial-correlation inference is invariant to the slab scale. Applies to GGM and mixed MRF models.
• bgm(delta =) sets the exponent of the determinant tilt on the continuous precision block. It defaults to NULL, resolved at fit time to 0.5 * log(p) over the continuous variables.
• bgmCompare(difference_family =) chooses the family of the prior on the group differences, "Normal" (the default) or "Cauchy" (0.1.6.3's fixed behaviour). It governs the pairwise-interaction differences and the main-effect threshold differences alike, and under difference_selection = TRUE it is the slab of the spike-and-slab.

**Inclusion inference.**

• Rao-Blackwellized inclusion probabilities and Bayes factors, through extract_posterior_inclusion_probabilities(estimator = "rb") and extract_inclusion_bf(). For each indicator update the sampler records the one-step draw J = gamma + (1 - 2 gamma) alpha (pre-move state gamma, birth/death acceptance probability alpha), stored per chain and per edge in fit$raw_samples$rb_inclusion alongside the raw indicator draws. Averaging J is a lower-variance estimator of the inclusion probability than the raw indic