select.explore() gains a
method = "BMA" option: Bayesian model averaging as
an alternative to the hard Bayes factor threshold
(method = "BF_cut"). For each edge the posterior is a
spike-and-slab mixture with mass P(H0 | Y) at zero and the
posterior under the alternative otherwise; the reported network weights
(pcor_mat_zero) are the median of this
mixture, computed exactly under a normal
approximation of the posterior of the Fisher-z partial correlation
(truncated at 0 for one-sided hypotheses). The result is therefore
deterministic and does not require stored posterior draws. The new
prior.prob.H0 argument (default 0.5) sets the prior
probability of the null. All alternatives are supported. For
alternative = "exhaustive" the mixture has three states – a
spike at zero (H0), a positive slab (H+), and a negative slab (H-) –
mixed by the posterior hypothesis probabilities;
pos_mat/neg_mat/null_mat then
classify each edge by the sign of the model-averaged median (so an edge
is null when the median is 0, which can happen even when H0 is not the
most probable hypothesis). Thanks to Joris Mulder for the
suggestion.select.explore() returns edge inclusion
probabilities (incl_prob): a matrix of posterior
edge inclusion probabilities q * BF / (q * BF + 1 - q) with
prior inclusion probability q = 1 - prior.prob.H0
(BF_10 for "two.sided", BF_20 for
"greater"/"less", and
1 - P(H0 | Y) for "exhaustive").prior.prob.H0 now affects the selected graph
under method = "BF_cut", not only under
method = "BMA". An edge is selected when its posterior
inclusion probability exceeds BF_cut / (BF_cut + 1), and is
called null when the posterior probability of the null hypothesis
exceeds that cutoff. With the default prior.prob.H0 = 0.5
this is identical to the previous Bayes factor rule
(BF > BF_cut) for "two.sided",
"greater" and "less". For
alternative = "exhaustive" the null keeps
prior.prob.H0 and the two directional hypotheses split the
remainder equally (0.5 / 0.25 / 0.25 by default); because
BF_1u + BF_2u = 2 the inclusion probability then equals the
one of "two.sided", so both alternatives select the same
edges, with "exhaustive" labelling each selected edge
positive or negative by the larger of the two directional probabilities.
alternative = "exhaustive" with
method = "BF_cut" now also returns
pcor_mat_zero. Previously "exhaustive" used
fixed equal 1/3 priors and thresholded each hypothesis separately, which
made BF_cut = 3 correspond to a two-sided Bayes factor of
1.5.truncnorm added to
Imports: required for the truncated-normal draws
used in the one-sided BMA alternatives.explore() no longer samples the prior by
default: the Bayes factors only need the prior sd of the
Fisher-z partial correlations, which is now computed analytically
(numerical integration of the marginal prior
rho ~ 2 * Beta(delta/2, delta/2) - 1, independent of
p) and stored as prior_sd_z. Draws from the
joint prior (prior_samp) are returned only with the new
argument store_prior_draws = TRUE (default
FALSE), which halves time and memory use for large networks
by default. select.explore(),
ggm_compare_explore() and bggm_missing() use
prior_sd_z; older explore objects still work
in select(). Bayes factors change slightly, because the
previously sampled prior sd was affected by the matrix-F approximation
(e.g. about 0.69-0.70 instead of 0.684 for
prior_sd = 0.5).explore() gains store_post_draws
(default TRUE): all samplers used by
explore() now keep running posterior summaries of the
partial correlations (post_samp$pcor_mat,
pcor_sd, z_mean, z_sd). With
store_post_draws = FALSE the p x p x iter
arrays post_samp$pcors and post_samp$fisher_z
are not stored, so memory use no longer grows with the number of
iterations. select() (all methods and alternatives) and
summary() work from these summaries when the draws are not
stored; functions that need the draws (posterior_samples(),
convergence(), coef(),
pcor_to_cor(), predict(),
predictability(), posterior_predict(),
constrained_posterior()) stop with an informative error,
and bggm_missing() always stores the draws.explore() gains burnin (default
50) and thin (default 1): after
burnin iterations, iter post-burn-in
iterations are run and every thin-th draw is stored; the
running posterior summaries use all post-burn-in iterations.
explore() objects no longer store the burn-in draws:
post_samp$pcors, fisher_z, beta
and thresh now contain only the kept draws (previously
iter + 50, including 50 burn-in draws). Downstream
functions locate the post-burn-in draws with an internal helper, so
objects created with earlier versions keep working. Objects from
estimate(), confirm() and the other functions
are unchanged.explore(iter = ) is now the number of
post-burn-in ITERATIONS, not the number of draws kept after
thinning. Previously explore(iter = i, thin = t) ran
burnin + i * t iterations and stored i draws,
so raising thin silently multiplied the run time by
t; it now runs burnin + i iterations and
stores ceiling(i / thin) of them, returned in the new
n_draws element. iter in the fitted object is
the iteration count and n_draws the stored-draw count. With
the default thin = 1 nothing changes.seed = NULL no longer re-seeds the
RNG: explore(), confirm(),
ggm_compare_confirm(), ggm_search() and
var_estimate() called set.seed(seed)
unconditionally before the guarded
if (!is.null(seed)) set.seed(seed). With the default
seed = NULL the first call was set.seed(NULL),
which re-initialises the RNG from the current time and process ID,
discarding any seed the user had set and making repeated calls on the
same data irreproducible even inside set.seed(). The
unguarded call has been removed, so with seed = NULL these
functions now use the ambient RNG stream and are reproducible under
set.seed(), as estimate() already was. Results
of existing scripts that relied on the implicit re-seeding will
change.alternative = "exhaustive", the
positive/negative Bayes factors were referenced to the null model
H0 instead of the unrestricted model Hu, so
each carried an extra factor of the two-sided Bayes factor and
double-counted the two-sided evidence. All three Bayes factors are now
referenced to Hu per Eq. 9 of Williams & Mulder (2019).
Affects post_prob and the derived
null_mat/pos_mat/neg_mat under
both method = "BF_cut" and
method = "BMA"."greater",
"less" and "exhaustive" selection branches of
select.explore(): the "greater"
branch averaged a hardcoded 3x3 prior mask
(upper.tri(diag(3))), giving wrong Bayes factors whenever
the number of variables was not 3, and the
"less"/"exhaustive" branches averaged the
(near-zero) matrix diagonal into the prior standard deviation. All
branches now use the off-diagonal (edge) prior sd, matching the
"two.sided" branch.summary() for
select.explore() with
alternative = "less": the "less"
summary branch produced an empty Relation column
(mat_names[upper.tri(mat_names)] on an already-flattened
vector), causing a “differing number of rows” error. The greater/less
summaries now share one correct branch.select.explore() now uses all post-burn-in
draws for the posterior mean/sd (previously the last 50 stored
draws were dropped), and no longer returns NA on the
diagonal of the exhaustive
null_mat/pos_mat/neg_mat (where
the posterior sd is 0).bggm_missing() dropping the wrong column
with mice >= 3.17.0 (#2):
bggm_missing() removed the .id column from
mice::complete(action = "long") by position (column 2),
which was correct only before mice 3.17.0. Since
mice 3.17.0 places .imp/.id in
the last two columns, this stripped a real data column and let
.id leak into the model, corrupting the data passed to
estimate()/ explore(). The column is now
removed by name, which is robust to mice’s column order.
Reported by @LilyTeesson.bggm_missing() now pools the posterior draws of
the imputed data sets correctly: it stacked all draws of each
fit, including their 50 burn-in draws, and set iter to
iter * m + 50, so the methods that use draws
51:(iter + 50) included the burn-in draws of imputations 2
to m and indexed beyond the stored draws; summaries came from the first
imputation only. The pooled object now contains the 50 burn-in draws of
the first fit followed by the post-burn-in draws of all fits, with
iter = iter * m. beta and thresh
are pooled whenever present, and pcor_mat (and, for
explore, the posterior summaries) are recomputed from the
pooled draws. Also works for m = 1.explore() for large or n < p
networks: the starting value is regularized
(solve(cov(Y) + 0.1 I)), and the matrix-F prior now uses
epsilon = min(0.01, 1 / (10 p)) instead of a fixed 0.01, so
that nu = 1 / epsilon stays well above p - 1
(required for a proper prior, Williams & Mulder, 2020). Results
change slightly for p > 10.delta from
prior_sd).rref_ei not found error: Added
simple_rref() function in helpers.R to replace
commented-out pracma::rref() call. This fixes a crash in
create_matrices() when validating constraint matrices for
hypothesis testing with multiple groups.ggm_search() crash: The C++
search function now handles edge cases where the adjacency
matrix becomes all zeros or all ones, preventing “sample more elements
than in x” errors.qplot() deprecation warning:
Replaced deprecated qplot() with ggplot() +
geom_density() in plot_prior().ggm_search() proposal-set shadowing bug
and BIC creep: The C++ search() function had a
variable-shadowing bug where zeros and
nonzeros (pools of candidate edges) were re-declared
instead of re-assigned after accepting a move, freezing them at the
starting graph’s configuration for the entire run. This caused
systematic BIC drift under probabilistic acceptance and made reversing
bad moves impossible. ggm_search() now implements a proper
Metropolis-Hastings sampler with a birth-death Hastings correction as
the default (probabilistic = TRUE); the greedy
deterministic hill-climb is still available via
probabilistic = FALSE. The burn_in parameter
is restored and now actually applied to discard pre-burn-in samples
before computing the Bayesian Model Averaging (BMA) solution
(probabilistic search only). Note: the greedy hill-climb was found to be
nearly non-functional on realistic test cases, typically accepting only
a single edge flip out of thousands of attempts.explore() (also used by
estimate(), confirm() and the
ggm_compare_*() functions) and by
var_estimate(): Theta_mcmc,
cors_mcmc and Sigma_mcmc (each p x p x iter)
were allocated but never filled or returned in
Theta_continuous, sample_prior,
mv_continuous, mv_binary,
mv_ordinal_albert, copula and
var; copula also allocated an unused n x p x
iter array of latent data. This reduces memory use substantially for
large networks (e.g. mv_continuous, used with
formula, allocated five p x p x iter arrays and now
two).missing_copula (mixed data with
impute = TRUE) no longer returns
post_samp$Y_collect: this n x p x iter array was
never filled (all zeros) and was not used anywhere in BGGM.rlang to Imports for proper use of
.data pronoun in ggplot2 aesthetics.bggm_fast.cpp for clarity and
efficiency: Removed two functions unreachable from R
(mv_ordinal_cowles, trunc_mvn; ~390 lines).
Extracted cov_to_cor() and
conditional_normal_params() helper functions, eliminating
17 duplicated blocks and a matrix inversion that was being computed
twice per variable per MCMC iteration across 9 sampler functions.
Removed dead commented-out debug code. No behavior change; file reduced
from 3199 to ~2850 lines.ggm_search() test suite and added
regression tests: Corrected tests that were silently exercising
non-existent parameters (bma = TRUE instead of
bma_mean, a start = argument that doesn’t
exist, gamma = instead of prior_prob,
prior_sd = which isn’t a parameter at all) and a
plot.ggm_search() test that could never fail since no such
method exists. Rewrote with correct parameter names and assertions on
the returned object’s fields, and added regression tests tied directly
to the bugs above (probabilistic search must accept a healthy fraction
of proposals; bma_posterior() must not error on the new
default’s output). Added tests/testthat/test-search_cpp.R
exercising the C++
search()/bic_fast()/hft_algorithm()
functions directly, including a check for systematic BIC drift across a
long run.current_thresh, candidate_thresh,
thresh_mat, c_thresh_mat, and
thresh_mcmc now begin from the same baseline after
initialization to avoid drift in Metropolis proposals.Removed NPM library to avoid CRAN compiler errors
search Functionstart_adj (the
maximum likelihood solution) as the starting point to avoid inefficient
sampling.zeros and nonzeros is now
adapted to the newly accepted adjacency matrix
(adj_s) rather than the static start_adj.find_ids(adj_mat) instead of
find_ids(start_adj) to ensure edge modifications are
correctly tracked after acceptance.Example in bggm_missing.R reintroduced; Was removed due to irreproducible CRAN error.
Addressed CRAN check error when building vignettes (removed example in bggm_missing.R – will put it back once it’s accepted to CRAN).
Downgraded required R version to 4.0.0
conv_to<>::from with as_scalarprior_sd: Adjusted computation of delta. Also, changed
default value for estimation: sqrt(1/3) resulting in delta = 2. For
model testing default is more tight, at sigma_sd = 0.5,
resulting in delta = 3.prior_sd is now limited to range 0 – sqrt(1/2)BFpack dependency error fixed.
This version of BGGM included changes based on the JOSS reviews: see here for the overview and here for specific issues.
BGGM was almost completely rewritten for version
2.0.0. This was due to adding support for binary, ordinal,
and mixed data, which required that the methods be written in
c ++. Unfortunately, as a result, lots of code from version
1.0.0 is broken.
Full support for binary, ordinal, and mixed data. This is
implemented with the argument type
roll_your_own: compute custom network statistics
from a weighted adjacency matrix or a partial correlation
matrix
pcor_to_cor: convert the sampled partial correlation
matrices into correlation matrices.
zero_order_cors: compute zero order
correlations
convergence: acf and trace plots
posterior_samples: extract posterior
samples
regression_summary: summarize multivariate
regression
pcor_sum: Compute and compare partial correlation
sums
weighted_adj_mat: Extract the Weighted Adjacency
Matrix
pcor_mat: Extract the Partial Correlation
Matrix
Five additional data sets were added.
ggm_compare_ppc: added option for custom network
statistics
Added option to control for variables with
formula
A progress bar was added to many functions
Initial CRAN release