This release rebuilds most of the package on top of 0.1.6.3.
bgm() now fits Gaussian graphical models and mixed
discrete-continuous networks alongside ordinal ones; priors are supplied
as objects instead of loose scalars; inclusion inference is
Rao-Blackwellized throughout; bgmCompare() has been
reworked; and a set of checking tools (verdicts(),
calibration_check(),
prior_sensitivity_check(), network and edge plots) is new.
Several defaults changed, some arguments were removed, and pairwise
parameters are reported on a different scale, so a script written
against 0.1.6.3 can run and give different numbers. Those changes are
listed first.
Throughout this file and in the package’s output, “log” means the natural logarithm.
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 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 indicator average
(estimator = "raw"). The sampler also accumulates the
inclusion odds on the acceptance-probability scale, so
extract_inclusion_bf() stays finite down to log-acceptances
of about -745: edges whose raw average saturates at 0 or 1 because the
chain never flipped in a finite run still receive a finite, data-driven
Bayes factor. The prior inclusion odds are divided out (from
extract_prior_inclusion_probabilities() for
bgm(), from the exchangeable difference prior for
bgmCompare()), so the reported value is a Bayes factor
rather than posterior odds, and equals the posterior odds only at a
prior inclusion probability of 1/2. Because the Rao-Blackwellized draw
is continuous, the standard MCSE/ESS/split-R-hat machinery applies to
inclusion inference. Available for bgm() (ordinal,
continuous, and mixed) and bgmCompare() (difference
indicators). The estimator inherits the chain’s mixing and does not
rescue a chain that failed to explore the model space.
extract_inclusion_bf(log =) selects the scale. The
default log = FALSE returns the inclusion Bayes factor
itself, matching the function’s name; log = TRUE returns
its natural logarithm. The accumulators are exact on the log scale
everywhere, while the Bayes-factor scale saturates at double precision:
an entry whose log exceeds about 709.78 nats (a Bayes factor beyond
about 1.8e308) is reported as +Inf under
log = FALSE even though its log-scale value is finite, and
an entry whose log is -Inf becomes 0.
Workflows that must separate such extreme evidence should use
log = TRUE. The argument applies to both the
bgm() and the bgmCompare() method, including
the stochastic-block difference path, whose return remains posterior
odds.
extract_prior_inclusion_probabilities() returns
prior edge-inclusion probabilities in the same matrix layout as
extract_posterior_inclusion_probabilities(), for
prior/posterior inclusion-odds computations. Under the joint
spike-and-slab prior on a continuous block the graph marginal is
reweighted by the per-graph normalizer (positive-definite-cone mass
shaped by the determinant tilt), so continuous-continuous edges do not
keep the edge-prior marginal at any delta — with a uniform
Beta-Bernoulli hyperprior at 3 variables the prior edge probability is
about 0.37 rather than 0.5, and a fixed
bernoulli_prior(0.5) at delta = 0 lands near
0.27. The values are read from the cached correction table
(bernoulli_prior(), beta_bernoulli_prior()) or
estimated by a prior-only chain with the fit’s own correction settings
(sbm_prior(), cached on the fit). Mixed models report
per-edge-class values (discrete-discrete, continuous-continuous, cross);
ordinal models use the analytic edge-prior marginals, including the
exchangeable partition mixture for sbm_prior().
Checking a fitted model.
verdicts() reads each edge (or difference)
indicator’s inclusion Bayes factor as a three-way verdict — evidence of
presence, undecided, evidence of absence, at a threshold and its
reciprocal — and flags the verdicts that a rerun of the sampler could
change. The flag is the union of two standard errors of the logit
inclusion probability: the Jeffreys-smoothed two-state model of the
binary indicator chain, which stays defined when the indicator never
flips, and the delta-method transform of the Rao-Blackwellized Monte
Carlo standard error. An edge is fragile when either places a verdict
boundary within two standard errors of the estimated evidence. The rule
comes from a known-truth calibration study of 37,010 graded edge-fits
across ordinal, binary, and Gaussian graphical models, in which every
one of the 66 verdict errors sat within 0.58 of a threshold on the log
Bayes factor scale and no edge further out was ever misclassified: over
those fits the two-state standard error alone recalled 0.742 of the
errors and the Rao-Blackwellized one 0.939, while the union recalled all
of them at a 3.0% false-alarm rate and transferred across model types.
Classification goes through the same internal routine as
prior_sensitivity_check(), so the two classify identically.
Printing reports the verdict tally and, when any verdict is fragile,
says so. Every arm of that study fitted a single network, so the
operating point belongs to the edge indicators of bgm(); on
bgmCompare() difference indicators the flag still marks
verdicts near a boundary, but its error-catching rate and false-alarm
rate there are unmeasured, and the print method says so rather than
borrowing the single-network numbers.
calibration_check() fits the isotonic
(pool-adjacent-violators) reliability curve of the model’s conditional
predictions, one panel per variable: pooling across variables lets a
variable predicted too high and one predicted too low cancel, and the
pooled curve then tracks the diagonal while neither variable does. The
consistency band resamples each case’s category from its own predicted
distribution rather than resampling the cumulative threshold events
independently, because those events are nested within a case and
treating them as independent understates the band. Continuous variables
are covered through the probability integral transform:
u = F(y | rest) is uniform exactly when the conditional
predictive distribution is right, so the panel is the empirical
distribution function of the u values against the uniform
diagonal. F is the predictive mixture over
ndraws posterior draws rather than a plug-in Gaussian at
the posterior mean, so parameter uncertainty sits inside the
distribution the observation is transformed by; the two panel kinds
therefore differ in what they condition on, and the kind
column of both returned tables records which construction produced each
row. The continuous band is not resampled from the model, because under
the transform the null is Uniform(0, 1) whatever the
conditional density was: it is the simultaneous envelope of the
empirical distribution functions of n independent uniforms.
Both panel kinds live on the unit square with the diagonal as the
calibrated reference, so a mixed fit produces one figure and one summary
table. Evaluated on the fitted data by default; newdata
makes it out-of-sample. In-sample the observation also entered the
parameters it is judged against, which makes the transform mildly
under-dispersed and the check conservative. A bgmCompare()
method is available and reports per-group panels.
prior_sensitivity_check() recovers the continuous
inclusion-Bayes-factor curve of each edge across the interaction slab
scale, and classifies how every edge’s verdict behaves along it. The
curve is anchored at a handful of fixed-scale fits (log-spaced
anchors, default 0.4 / 0.63 / 1 / 1.6 / 2.5
times the chosen scale) and filled in between anchors by importance
reweighting: a fixed-scale fit is reweighted to a nearby scale with
per-draw slab-density ratios over the included edges, the likelihood
cancels, and no normalizing constant is involved. The 1x
anchor is the original fit itself and is never refit, so the
chosen-scale verdicts the check reports are exactly the analysis already
run. Each display point pools every anchor that clears
ess_floor (default 400) there, weighting each anchor’s
inclusion-probability estimate by its inverse variance on the
probability scale and transforming the pooled probability to the
Bayes-factor scale; points where no anchor clears the floor are reported
NA rather than extrapolated, and non-overlapping anchor
radii warn to add anchors. Exactness is kept off the pooled curve and on
the anchor fits themselves: the per-anchor verdict columns and the
chosen-scale quantities come from those fits’ own Rao-Blackwellized
statistics, not from a pooled row. An edge is called scale-sensitive
only if its verdict changes along the curve and its
Bayes-factor swing exceeds a run-to-run noise band (the 95th percentile
of the spread between one anchor refit and a repeat of it), so boundary
edges that merely flip between reruns of the same prior are reported as
indistinguishable-from-wobble rather than as moves; an edge
that is threshold-relevant in one of the two refits and saturated in the
other has a censored rather than an infinite spread and is left out of
that percentile, with the printed report counting the exclusions. Each
anchor passes a convergence gate (median continuous split-R-hat, bulk
Rao-Blackwellized inclusion R-hat, and E-BFMI for NUTS) before it enters
the curve. The check works on any bgm() fit with edge
selection; on Wenchuan it costs about one original fit with warm-started
NUTS refits (a refit_sampler argument, default
"same-as-fit", carries the adapted NUTS step size and mass
matrix).
prior_sensitivity_check(vary =) names which prior
the curve moves, and the printed report says which. On a continuous or
mixed fit the interaction slab scale and the prior on the precision
diagonal are tied through the standardized frame (raw rate
= eta / s), so a sweep of the slab has to hold one of the
two fixed, and the two choices answer different questions:
vary = "slab" holds the raw diagonal rate at the fitted
value and moves the interaction prior alone, while
vary = "slab-and-diagonal" holds eta fixed and
lets the raw rate follow, rescaling the prior’s overall size with its
shape held fixed. The default vary = "auto" follows the
frame the fit itself used — the joint sweep when the diagonal prior was
given as eta, the slab alone when it was given as a raw
rate. On a five-variable Gaussian graphical model the two
modes give visibly different curves, so the choice is not a formality.
Discrete fits have no precision diagonal; vary is accepted,
has no effect there, and the report omits the line.
prior_sensitivity_check() works on
bgmCompare() fits, tracing each difference indicator’s
inclusion Bayes factor across the difference slab scale.
bgmCompare() gives the pairwise and the main-effect
difference families one difference_scale, so a single curve
covers both. The anchored machinery, the convergence gate, the noise
band, and the reporting are the same as for bgm(); what
differs is that one difference indicator gates several parameters (a
pairwise difference in every group contrast, a main-effect difference
across that variable’s whole threshold block), so the importance weight
sums the difference-slab density ratio over the gated parameters while
the inclusion probability is still reported per indicator. Reweighting a
fit at one scale to another was checked against real refits at that
scale: on Wenchuan at a doubling of the scale the reweighted inclusion
probabilities land within 0.016 of a refit’s, against a 0.007 run-to-run
spread between two refits at that same scale, at an importance effective
sample size of 2536 out of 12000 draws. Indicators the sampler never
updated carry no verdict at any scale, and the printed report counts
them out rather than reporting them as undecided. A stochastic-block
difference prior has no single marginal inclusion probability and is
refused, since the curve would report posterior odds rather than Bayes
factors.
New vignette “Checking your fitted model”:
verdicts() and the verdict-encoded network, then
calibration_check(), then a worked example of building a
predictive display directly on simulate() — the joint-level
check is an analyst’s reading of a display, not a pass/fail statistic,
so the vignette shows how to roll one rather than calling a blessed
function. The worked statistic is the sum-score distribution against its
pointwise predictive band, chosen because it lies outside the model’s
sufficient set: pairwise dependence statistics sit close to that set, so
how they read depends on how the estimator is anchored rather than on
whether the model fits.
Plots and centrality.
plot() on a bgm() fit draws the
edge evidence plot: three panels on one shared layout,
holding the pairs the data support, the pairs the data rule out, and the
pairs the data cannot decide. A single network drawing has to make every
pair either an edge or a blank, and a blank cannot say which of those
two a missing edge is; because bgm() returns an inclusion
Bayes factor for every pair, that choice does not have to be made. Each
panel is titled with what it holds and how many pairs are in it, and
with the rule that put them there, stated as a Bayes factor rather than
its logarithm. Only the first panel is weighted: line width is the
posterior mean pairwise association and colour carries its sign, because
that is where the effect sizes are; the colour convention is documented
in ?plot.bgms. The other two are drawn at uniform width,
dashed and dotted, because for those pairs the classification is the
result. The layout is computed once from every pair, so a node sits in
the same place in all three panels, and an all-absence fit is a result
rather than an error: the absence panel fills and the presence panel is
empty. The threshold is the one verdicts() uses, so the
picture and the reporting table state the same thing. A fit run without
edge selection has no inclusion Bayes factor and nothing to split its
pairs by, so plot() draws one weighted panel titled
Edge weights — the title says which channel the figure is
drawn in, so a wide line is never ambiguous between a large association
and strong evidence for one. type = "centrality" draws the
centrality display instead. Layout comes from qgraph, which stays a
suggested package; without it the method errors and points to
verdicts().
plot() on a bgmCompare() fit draws the
same three panels for difference evidence, with the supported panel
weighted by the posterior mean difference. The display draws for any
number of groups: a pair has a single inclusion indicator shared across
all K - 1 contrasts, so the three-way split is as well
defined for three groups as for two; what a pair does not have beyond
two groups is one magnitude, so the supported panel is weighted only for
two groups and drawn at uniform width otherwise, where the
classification is the whole of what it reports. A fit with no
differences left is a result rather than an error — “the groups do not
differ anywhere” is a common and correct finding. Main-effect
differences are not edges, and when
main_difference_selection = TRUE gives them their own
indicators their evidence rides on the nodes: each node wears a ring
filled to that difference indicator’s posterior inclusion probability —
a full ring is 1, half a ring 0.5. Under the default
main_difference_selection = FALSE those indicators are
never updated and no ring is drawn. type = "groups" draws
each group’s own posterior mean network on the same layout and pages
them, so a fit with more groups than max_panels is drawn a
page at a time rather than squeezed into one row; those panels use the
same colour-vision-safe sign pair as the rest of the package rather than
qgraph’s green/red default. With
difference_selection = FALSE, weights are the display: two
groups get one network of every pair, width the posterior mean
difference and colour its sign, titled Difference weights;
beyond two groups the weights of a pair are K - 1 numbers
whose honest weighted picture is the groups themselves, so
plot() draws the type = "groups" panels
without being asked for it. Read the magnitudes per group with
plot(fit, type = "groups") or
extract_group_params(). type = "centrality"
draws the centrality display, with group passed
through.
plot_edge_posterior() draws one edge the way JASP
draws a parameter: the posterior of the edge weight against the prior it
was updated from, on a numbered density axis. Under edge selection the
posterior drawn is the weight given the edge is included — a
real density integrating to one, which is what puts numbers on the y
axis — and the mass at zero moves off the plot region onto a probability
wheel filled to the posterior inclusion probability, matching the ring
encoding the bgmCompare() panels use for their main-effect
differences. The prior curve is computed in closed form from the fit’s
own interaction_prior rather than assumed, so a
cauchy_prior(scale = 2.5) fit and a fit at the default
normal_prior(scale = 1) do not look alike. The median, the
95% credible interval, and the natural-log inclusion Bayes factor are
printed as annotations in JASP’s positions.
What the panel shows follows the estimator the fit licenses. With
edge selection, the evidence is the Rao-Blackwellized indicator Bayes
factor, which is not a ratio of densities at zero, so no Savage-Dickey
ordinates are drawn — dots there would assert an estimator the package
does not use. Without edge selection there is no indicator, the
posterior of the weight is continuous, and Savage-Dickey is the
licensed estimator: the panel then draws the JASP figure exactly, with
both grey ordinate dots at zero, the ratio printed as a log Bayes
factor, and the wheel filled to BF/(1 + BF); what the two
shares are is documented in ?plot_edge_posterior rather
than tagged on the figure. The posterior ordinate is estimated by
stats::density(bw = "SJ") at zero; JASP’s own
implementations use a logspline fit, which is not adopted here because
it would add a dependency for one number. On a continuous block the
drawn prior is the slab the sampler evaluates, which is the marginal
prior of a precision entry only up to the positive-definite restriction;
?plot_edge_posterior says so. An edge that the data rule
out gets a figure rather than an error: the prior, the wheel (nearly all
pale) and the evidence are drawn. plot_edge_posterior()
refuses a bgmCompare() fit rather than drawing its baseline
pairwise effect as though it were one network’s edge; use
plot(fit) and verdicts(fit) there.
extract_centrality() evaluates a node centrality on
every posterior draw of the network and returns a draws-by-nodes matrix,
so the posterior distribution of a node’s centrality comes out of the
same model-averaged posterior as the edge weights, with structural
uncertainty included: an edge excluded at a given iteration contributes
exactly zero there. summary() reports the posterior mean,
credible interval, and the posterior probability of being the most
central node; plot() draws the ordered interval display.
measure = "strength" is the one accepted value for
now.
extract_centrality() works on
bgmCompare() fits. group = 1 gives that
group’s centrality, one row per posterior draw;
group = c(1, 2) gives the difference between two groups’
centralities, again per draw, so the credible interval is the interval
of the difference and answers whether the groups differ in a node’s
centrality directly, where two separately drawn intervals do not. Each
group’s network is rebuilt on every draw as
baseline + (P %*% differences) with the fit’s own contrast
projection, not from posterior means, which is what carries the
uncertainty through; averaged over draws it reproduces
extract_group_params() to machine precision. A draw in
which a difference indicator is excluded gives both groups the same edge
weight and contributes exactly zero, so a difference carries a point
mass at zero and summary() reports p_positive
(the probability the first group’s centrality is higher) rather than
p_most_central, with the remaining mass on no difference at
all. plot() marks zero.
New extractors.
extract_main_effects(): main effect samples — category
thresholds, continuous means, and the precision diagonal — replacing
extract_category_thresholds().extract_precision(): posterior precision matrix samples
from GGM and mixed models.extract_partial_correlations(): posterior partial
correlation samples from GGM and mixed models.extract_log_odds(): log-odds for discrete pairwise
interactions.Prior-only sampling.
sample_graph_prior() draws edge-inclusion
indicators, together with any edge-prior hyperparameters, from the graph
level of the spike-and-slab prior. Under the hierarchical specification
the draw is ancestral and exact; under the joint specification it runs
the zero-data prior chain, so the draws carry the per-graph normalizer
tilt. Hyperparameters can be fixed instead of sampled
(theta for Bernoulli and Beta-Bernoulli priors,
allocations plus block_probs for the
Stochastic-Block prior).
sample_sbm_prior() gives ancestral draws of block
allocations and pair-inclusion probabilities from the MFM-SBM edge-prior
hyperprior.
sample_ggm_prior() runs a prior-only GGM chain. It
accepts update_method = "gibbs" for the
spec = "joint" prior chain, using the conjugate row and
edge updates instead of adaptive Metropolis; edge_prior
objects (bernoulli_prior(),
beta_bernoulli_prior(), sbm_prior()),
returning sampled allocations under the block model; and
spec = "hierarchical".
Edge priors on continuous and mixed models.
Corrected inclusion-probability updates for
beta_bernoulli_prior() on continuous (GGM) models: under
the determinant-tilted precision prior, the conjugate Beta update omits
a normalizing-constant factor and biases the sampled inclusion
probability toward sparsity (at 5 variables with a uniform hyperprior
its prior mean lands near 0.37 instead of 0.5). The update now draws
from the corrected conditional using a table built from the prior
distribution at the first fit of a model configuration and cached on
disk (tools::R_user_dir("bgms", "cache"); a one-time cost
of the order of minutes, announced when verbose = TRUE).
The sampled inclusion probability is returned per chain in
fit$inclusion_parameter_samples.
Corrected stochastic block model updates for
sbm_prior() on continuous (GGM) models, using the same
cached table: the block-probability draws, the block-allocation weights,
and the new-cluster weight all carry the normalizing-constant
correction. Without it the sampled partition collapses toward one block
(at 20 variables the prior mean number of blocks lands near 1.0 instead
of 1.87); with it a prior-only chain reproduces the model’s partition
prior (total variation 0.01-0.02 at 5-20 variables).
The normalizing-constant correction extends to mixed models with
beta_bernoulli_prior() or sbm_prior(): the
determinant tilt acts on the continuous precision block, so the
correction table is built for the continuous variables and the
block-structure corrections read continuous-continuous edges only. With
fewer than two continuous variables no edge is tilted and the plain
conjugate updates apply unchanged; with exactly two, the single tilted
pair supports the inclusion-probability correction but not the
block-model slope curve, so sbm_prior() warns and keeps the
plain conjugate updates there. Prior-only mixed chains reproduce the
Beta hyperprior on the inclusion probability and the partition prior on
the number of blocks.
The correction table announces itself only when it actually
builds (a cache hit is silent) and states that the build is one-time and
cached. In interactive sessions it draws a progress bar for both serial
and parallel builds. The bar follows display_progress (the
same control as the sampler’s bar), not the advisory
bgms.verbose flag.
Under precision_graph_prior = "joint", a fit with
edge selection on a continuous precision block reports once that the
realized edge-inclusion prior is pi(Gamma) * Z(Gamma), the
edge prior reweighted by the per-graph normalizer, rather than the
nominal edge prior. This holds for every edge prior: at three variables
a uniform beta_bernoulli_prior(1, 1) realizes about 0.37
and a fixed bernoulli_prior(0.5) at delta = 0
realizes about 0.27. With a learned inclusion probability
(beta_bernoulli_prior(), sbm_prior()) the
corrected hyperparameter update keeps the update coherent with the joint
model but does not restore the nominal prior; with a fixed inclusion
probability nothing absorbs the tilt. The message follows
verbose (and the advisory bgms.verbose flag)
and points to extract_prior_inclusion_probabilities() for
the realized prior and to
precision_graph_prior = "hierarchical" for the nominal
one.
The hierarchical precision-graph prior (continuous and mixed models).
precision_graph_prior = "hierarchical" is the
default on continuous and mixed models. It composes the edge prior and
the precision prior as p(Gamma) p(K | Gamma) with
p(K | Gamma) normalized per graph, so the graph marginal is
exactly the edge prior; under
precision_graph_prior = "joint" it is that prior reweighted
by the per-graph normalizer. Both specifications are fully supported;
the default is the one whose graph marginal is the prior you wrote down.
Each edge move evaluates the normalizer ratio with a fast per-component
surface approximation, built once per analysis at the fixed
(eta, delta) from block-Gibbs anchors. It requires
edge_selection = TRUE and a normal_prior() or
cauchy_prior() interaction prior; the precision scale prior
may be exponential_prior() or a gamma_prior()
of any shape. On mixed models the specification normalizes the
continuous block p(K_yy | Gamma_yy): the approximation
enters the continuous-continuous edge moves only, and at least two
continuous variables are required. The same specification is available
in sample_ggm_prior(spec = "hierarchical"). In simulation
checks with data (50 variables, n = 25, Beta-Bernoulli(2, 4)) the
posterior inclusion probabilities are unbiased to within +0.002. Chains
sampled from the prior alone, with a sampled inclusion probability and
many variables, can fall outside the approximation’s validated range;
the trust gauge flags affected chains.
The argument is accepted where the choice between the two
specifications is vacuous rather than erroring. The specifications
differ only in how p(K | Gamma) is normalized across
graphs, so they differ only where the sampler moves between graphs. With
edge_selection = FALSE the graph is fixed and the two
coincide exactly; the argument is accepted silently, and the correction
surface and the trust gauge are skipped rather than built and left
unused. With no continuous precision block — an ordinal model, or mixed
data with fewer than two continuous variables — there is no
K for the argument to refer to; the fit is accepted, and a
message reports this when the argument was named and
verbose = TRUE. Inheriting the default is silent: since
"hierarchical" is the default, the value reaches every fit,
and the message reports a user’s choice rather than a default nobody
made. Both cases run the joint path, which applies no correction on
exactly these configurations either, so the posterior is identical to a
"joint" fit from the same seed. A
beta_prime_prior() interaction slab is still rejected where
the choice is meaningful — a continuous precision block under edge
selection — and is tolerated where the argument has no
referent.
The per-edge normalizer ratio is corrected by a theta-independent
absolute-moment surface. The surface predicts the per-component first
and second spectral moments for both mediating-block families
(common-neighbour clusters and bipartite bridges), and is deployed by
decomposing each block into disjoint components. Its constants are built
at unit slab scale, where the normalizer ratio is invariant to the slab
scale, so results do not depend on the user’s scale choice, and a
cauchy_prior() slab gets its own constants rather than
reusing the Normal ones. Validated against a per-component gold
reference across a graph-density sweep: about 0.003 nats to gold,
against about 0.036 for the additive kernel it replaces, up to about
1000x tighter in the dense regime. The surface deploys at every
standardized rate eta, with gold gates recorded at
eta = 1 and eta = 2.
Where the surface deploys, and what serves the cells it does not:
gamma_prior() diagonal shapes from 0.5 to 10 are served
by the surface. It is scored against a block-Gibbs reference at shapes
0.5, 1, 2, 3 and 5 — every interior cell inside the 0.003-nat envelope,
against 0.016 to 0.40 nats for the additive kernel — and the interior of
that range is interpolated, not measured. Anchors at a non-unit shape
are run at a raised sweep budget, resolved by matching the measured
anchor spread to the shape-1 reference.eta) of 1 and 2;
past 2 the same route deploys and a note says plainly that the bound
does not cover the cell.delta across 6 to 120 variables, the boundary runs from
about 12 to 32 common-neighbour variables; the collapse is out of reach
below roughly 16 variables and becomes reachable from about 20 to 30
upward. A fit on which it happens says so, driven by the counters below,
and reports the largest common-neighbour block involved.
?summarize_zratio_gauge documents the mechanism and the
measurement.The anchored size range reaches 80 variables, so the reachable
mediating block of an ordinary large fit sits inside it rather than in
the extension zone; an anchor tier is placed at the size cap itself
whenever the surviving top tier is more than two sizes below it. Past
the anchored range the surface continues along its own boundary slope in
log-size rather than freezing at the range’s edge. Scored against a
block-Gibbs reference on single-component blocks of 90 to 150 variables
against an 80-variable anchored range, the extension holds a median of
0.0006 nats and at most 0.0011 for common-neighbour blocks, and a median
of 0.0043 and at most 0.0060 for bipartite ones, where freezing reaches
0.060 and 0.095 and the additive fallback reaches 0.21 and 0.31. Inside
the anchored range the surface is at 0.003 or better, so the extension
is bounded and much better than the alternatives rather than validated
to the in-range envelope, and the bipartite bound in particular sits
about twice outside it. The absolute moments grow with component size,
so a fitted slope pointing downward at the range edge is a fit
pathology; it is floored at zero and counted in the per-chain
n_slope_floor diagnostic.
The certified range of the tabulated constants is [0.5, 20]; the fixed quadrature grids were scored at shapes 2, 10, 12, 15 and 20 at standardized rates 1 and 2 against references sharing no quadrature with them, worst deviation 8.6e-08 against a 1e-07 to 1e-06 working tolerance. A shape above 20 warns.
The surface build is cached (session and disk, keyed on the fit
cell and size cap, following
options(bgms.correction_table_cache)), so repeated fits of
the same configuration reuse it. A first build uses the fit’s own
cores (the chains have not launched yet; the result is
identical for any core count) and can be overridden with
options(bgms.zratio_surface_cores); on Windows a socket
cluster stands in for forking on large builds. It announces itself when
it runs, reporting the size cap and the number of cores it uses, and
reports the elapsed time when it finishes — this accounts for the pause
before the chains start on a first hierarchical fit. A cache hit is
silent, and the announcement follows the advisory
bgms.verbose flag. The build costs about 24 s serial and 7
s on four cores at the full 80-variable range at the exponential shape
(43 to 47 s serial at a non-unit shape); a fit on fewer variables sizes
the build to its own variable count.
A trust gauge runs by default on the deployed hierarchical path —
which is to say on a default continuous or mixed fit with edge selection
— with options(bgms.zratio_gauge_sweeps = 0L) as the off
switch; sample_ggm_prior() keeps its own
zratio_diagnostics argument. While the sampler runs, the
gauge redoes a subset of each chain’s edge add/remove decisions with the
exact calculation and records flip_rate, the fraction of
decisions that would have come out differently; a chain whose flip rate
exceeds 1% is flagged. A second channel, harm_pred,
projects the measured approximation error onto the inclusion-probability
scale using the chain’s own edge sensitivities and the edge-prior
feedback amplification, and flags when the projected distortion exceeds
0.02 — this catches a consistent error that shifts the recovered network
without changing individual edge decisions, which flip_rate
cannot see at chains whose decisions are far from their accept/reject
boundaries. It is computed for Bernoulli and Beta-Bernoulli edge priors.
The gauge is silent when clean and prints only flagged chains; the
per-chain summary is attached as fit$zratio_diag (and as
$zratio_diagnostics on sample_ggm_prior()
output). Its cost is fixed per chain rather than proportional to
iter: two assessment sweeps audit a capped number of edge
moves, which is nothing on a sparse posterior and about 5-10 seconds per
chain on a dense posterior at 100-200 variables. That is negligible
against a production-length fit and noticeable on a short exploratory
one, which is what the off switch is for. See
summarize_zratio_gauge() and the diagnostics vignette,
which documents the audit as a sample: it resolves coherent error, which
is what harm_pred projects and what shifts a recovered
network, and does not resolve rare edge-specific failures.
The advice printed for a flagged chain is a ladder. It reports
the measured error share, the mediating-block sizes the audit covered
(block_lo, block_hi) and how many of the
chain’s non-trivial moves were audited; then advises raising
options(bgms.zratio_gauge_sweeps) and refitting to resolve
whether the signal is real; and only then mentions
precision_graph_prior = "joint", stating that it targets a
different model, whose graph marginal is the edge prior reweighted by
the per-graph normalizer rather than the edge prior itself. The gauge
never switches specification on its own.
Per-chain diagnostics in fit$zratio_diag carry
flip_rate, amplification,
harm_pred, harm_flag, the
se_mcse/se_se uncertainty of the error
estimate, the audited block sizes (block_lo,
block_hi), the retained-sweep extrapolation counters
(n_pred_retained, n_extrap_retained,
max_extrap_size_retained), n_slope_floor,
n_isolated (edge evaluations taking the isolated-edge route
above shape 10) and n_collapsed /
max_collapse_size (the additive kernel’s zero-collapse,
counted at the moment-sign condition itself and excluding trust-gauge
sweeps). Post-fit notes are driven by those counters, so they report
what the chains did rather than what the specification intended, and
they reach a fit run with verbose = FALSE. The
extrapolation notice separates warmup from retained sweeps — the sampler
initializes from a complete graph, which puts every mediating block at
roughly the number of variables for the first sweeps — and says plainly
when no retained sweep extrapolated at all; trust-gauge sweeps run after
sampling and enter neither tally.
The “Getting Started with bgms” vignette has been rewritten
against this release. It opens on what the package is for – Bayesian
analyses of graphical models, all of them Markov random fields, with
variable_type selecting the model family – and on what an
analysis returns: a posterior over graphs and an inclusion Bayes factor
per pair, read as evidence of presence, evidence of absence, or
undecided. It states the current defaults
(normal_prior(scale = 1) for the interaction slab,
bernoulli_prior(0.5) for the edge prior,
"hierarchical" for precision_graph_prior) and
works through summary(), verdicts() on the
natural log Bayes-factor scale, and plot()’s three-panel
edge evidence display on nine Wenchuan items. It gained the bibliography
header the other vignettes carry, a bgmCompare() section,
and pointers to the package website and to easybgm. The
hand-thresholded median-probability graph drawn with
qgraph, the coef() walkthrough, and the
“What’s new in 0.2.0” section were dropped.
The “Model Comparison with bgmCompare” vignette now draws the
groups with plot(fit, type = "groups") – each group’s own
graph on one shared layout – in place of the hand-assembled
qgraph() call it built from
coef(fit)$pairwise_effects_groups, which drew one group
only. It gained the bibliography header the other vignettes carry, and
its introduction now states the inclusion Bayes factor and the three
verdicts verdicts() returns.
NUTS uses Stan’s multinomial candidate weighting (log-sum-exp of
H0 - h per leaf, biased progressive sampling at the top
level) in place of the Hoffman-Gelman slice variable. The two schemes
target the same posterior; the multinomial variant produces
lower-variance candidate selection and has been Stan’s default since
2017. User-facing output is unchanged apart from the new
accept_prob diagnostic.
NUTS Stage-2 warmup windowing matches Stan’s
windowed_adaptation::compute_next_window: when the window
after the next would overshoot the Stage-3a boundary, the current next
window is stretched to absorb the remaining Stage-2 budget instead of
emitting a small trailing window. This eliminates a disruptive
mass-matrix update and step-size reinit at the end of warmup and
improves dual-averaging convergence.
NUTS diagnostics include the per-iteration mean Metropolis
acceptance probability (fit$nuts_diag$accept_prob,
paralleling Stan’s accept_stat__) and a per-chain
mean_accept_prob summary. The diagnostics summary also
prints the warmup_incomplete flag (energy not stationary)
it already computed.
bgm(progress_callback =) and
bgmCompare(progress_callback =) take a function called as
sampling proceeds, so an embedding application (JASP) can drive its own
progress display.
Dropped coda from Imports; ESS and R-hat are
computed in C++ with on-demand (lazy) evaluation, replacing the eager
R-based computation of 0.1.6.3. $ and [[
accessors on fitted objects trigger that computation on first
access.
Refactored the C++ backend: unified model hierarchy
(BaseModel -> GGMModel /
OMRFModel / MixedMRFModel), shared NUTS
infrastructure, and fused log-posterior and gradient
computation.
On Windows with RcppParallel >= 6.0.0 (which
bundles oneTBB 2022), a multi-core run and a single-core run started
from the same seed are no longer bit-for-bit identical: the
second and later chains follow different sampling trajectories under the
new scheduler. Each chain still samples the same posterior — posterior
summaries agree to within Monte Carlo error — so results are
statistically equivalent, only not bitwise reproducible across core
counts. Other platforms and earlier RcppParallel are
unaffected.
extract_arguments() now reports
main_effect_indices for bgmCompare() fits: the
per-variable column blocks of the baseline main-effect parameters. The
layout was computed during fitting but kept internal, so there was no
supported way to map main-effect parameters back to variables.
Fixed the handling of ordinal categories in
bgmCompare(), a defect present in 0.1.6.3. The comparison
kept only the categories a variable was observed in every
group, and silently merged the rest into their neighbours: with two
groups observing categories 0-2 and 1-3 of a four-category item, both
groups’ data were folded onto two categories, however many observations
each category carried. Because groups being compared routinely differ in
the categories they use, this destroyed well-observed categories rather
than unused ones — a five-level variable could be reduced to a binary
one — and the affected variable’s pairwise parameters were then badly
overestimated, which also inflated the parameters of the variables it
connects to. bgmCompare() now keeps every category that any
group observes and merges none of them; only a category value that no
group observes at all is dropped, with the remaining categories
renumbered contiguously and a message() saying so. When a
kept category is empty in some group, that group’s threshold for it
rests on the prior rather than on data, so bgmCompare()
warns and names each variable, category, and group affected, and records
the per-group category counts in the fit
(extract_arguments(fit)$category_support). Blume-Capel
variables were never affected and are still exempt, because their
parameters are functions of the numeric category score. Any
earlier group comparison in which the groups did not observe exactly the
same categories should be refit; the pairwise estimates it
produced may be substantially too large. Comparisons whose groups shared
the same observed categories are unchanged.
Fixed the number-of-blocks summary for stochastic-block edge
priors, a mismatch a 0.1.6.3 user would have seen. The conditional p(K |
t) behind posterior_num_blocks placed a zero-truncated
Poisson prior on the number of components, while the sampler’s own
partition coefficients use the shifted Poisson (K - 1 ~
Poisson(lambda)). The mismatch reweighted the reported distribution by a
factor lambda/K across K, so the reported number of blocks was
systematically off from the one the sampler drew under. The summary now
uses the sampler’s convention and is unit-tested against the generative
partition prior.
Fixed a compilation failure on Alpine/musl that also affected
0.1.6.3: mrf_simulation.cpp used
tbb::global_control while relying on a transitive include
for <tbb/global_control.h>, which is not available on
all platforms. The header is now included directly.
Fixed a target_accept that was partly ignored, in
0.1.6.3 as well. The value reached dual averaging, but the step-size
heuristic that reruns after a mass matrix update used a hard-coded
0.625 instead, so a fit set to a different
target_accept was restarted at the wrong target every time
the metric was re-estimated. The heuristic now uses the requested
target.
Fixed NUTS acceptance-probability accumulation, present in 0.1.6.3 too: the top-level trajectory loop overwrote the Metropolis contribution with the last subtree’s value instead of summing across the full trajectory, so the signal driving dual-averaging step-size adaptation was biased and the adapted step size with it.
predict() now recodes newdata through
the category map stored with the fit rather than by subtracting each
column’s minimum. Under the old rule a sparse category coding — one with
gaps, such as values 1, 2, 4, 5 — was miscoded silently: the fit
collapses those values to categories 0, 1, 2, 3, but the minimum shift
mapped them to 0, 1, 3, 4, so every value above a gap was attributed to
the wrong category and the largest value fell outside the fitted
category range entirely, with no error and no warning.
simulate() correspondingly returns ordinal data on the
original category scale, so the simulate() to
predict() round trip stays on the scale the model was
fitted on.
bgm(standardize =) and
bgmCompare(standardize =) — see Breaking changes;
FALSE warns, TRUE errors.bgm() and
bgmCompare() — pairwise_scale,
main_alpha, main_beta,
inclusion_probability, beta_bernoulli_*,
dirichlet_alpha, lambda — in favour of the
prior objects interaction_prior,
threshold_prior, means_prior, and
edge_prior. Passing one warns and is translated.bgmCompare(difference_probability =) in favour of
difference_prior = bernoulli_prior(); the argument keeps no
default.extract_category_thresholds() in favour of
extract_main_effects().extract_ess(estimator = "mixt"), the transition-based
effective sample size of the binary indicator chain: it warns and
recomputes the quantity from the raw indicator draws, and the
n_eff_mixt column is gone from the fit summary’s inclusion
table. Three grounds. It carries no decision value: in a known-truth
calibration study its best error-catching threshold rejected 41% of
correct edge verdicts to catch 85% of the wrong ones, while a rule built
on the inclusion standard error caught all of them at a 2.3% cost. Its
own model is violated: the conversion from flip counts to an effective
sample size assumes a two-state first-order Markov chain, and the study
rejects first-orderness for 41.3% of edge-fits on transitions out of the
included state (3.2% out of the excluded state), asymmetrically. And it
is redundant: the precision of the inclusion probability and its Bayes
factor is carried by the Rao-Blackwellized n_eff and
mcse, while the switching frequency mirrors the inclusion
probability, so a low n_eff_mixt mostly restates that the
inclusion probability is extreme. The raw directional flip counts stay
in the table.plot_edge_posterior(binwidth =) is deprecated and
ignored. The panel expresses the weight as a density rather than as
probability per weight bin, so there is no bin to set a width for.extract_rhat(): extract R-hat convergence diagnostics
from fitted objectsextract_ess(): extract effective sample size estimates
from fitted objectsverbose argument: control informational messages; set
options(bgms.verbose = FALSE) to suppress globallysimulate_mrf(): standalone MRF simulation with
user-specified parameterssimulate.bgms(): generate observations from fitted
models (supports parallel processing)predict.bgms(): compute conditional probabilities P(X_j
| X_{-j})main_difference_selection argument in
bgmCompare(): control threshold difference selectionstandardize argument: scale Cauchy prior by response
score rangebaseline_category now stored in fitted object for
Blume-Capel simulation/predictionposterior_mean_indicator: now
correctly maps C++ row-major order to R matrices (#77)(c - ref) in pseudolikelihood denominatormrfSampler() → use simulate_mrf()bgm() and
bgmCompare() modelssummary(), print(), and
coef() methods for fitted objectsbgm() and bgmCompare() is
now reproducible when a seed argument is specifiedsummary() now integrates the functionality of the old
summary_SBM()main_difference_model = "collapse" settingbgmCompare() when
handling missing databgm(), the following arguments are deprecated:
interaction_scale → use
pairwise_scaleburnin → use warmupsave → no longer needed (all outputs are returned by
default)threshold_alpha, threshold_beta → use
main_alpha, main_betabgmCompare(), arguments related to difference models
are deprecated:
main_difference_model (removed without
replacement)reference_category → use
baseline_categorypairwise_difference_*, main_difference_* →
use unified difference_* argumentspairwise_beta_bernoulli_*,
main_beta_bernoulli_* → use unified
beta_bernoulli_* argumentsinteraction_scale → use
pairwise_scalethreshold_alpha, threshold_beta → use
main_alpha, main_betaburnin → use warmupsave → no longer neededextract_edge_indicators() → use
extract_indicators()extract_pairwise_thresholds() → use
extract_category_thresholds()$gamma (pre-0.1.4) and $indicator
(0.1.4–0.1.5) → replaced by $raw_samples$indicator$main_effects (pre-0.1.4) and
$posterior_mean_main (0.1.4–0.1.5) → replaced by
$raw_samples$main (raw samples) and
$posterior_summary_main (summaries)selection = TRUE, the burnin phase now runs
2 * burnin iterations instead of 1 * burnin.
This ensures the chain starts with well-calibrated parameter valuesThis is a minor release that adds some documentation and output bug fixes.
This is a minor release that adds some bug fixes.
This is a minor release adding some new features and fixing some minor bugs.
na.action option.edge_prior option.adaptive
option.interaction_prior option.threshold_alpha and
threshold_beta options.