BGGM

• **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.

Bug fixes

• **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.
• **Corrected the exhaustive posterior hypothesis probabilities**: for 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".
• **Corrected the prior sd in the "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.
• **Fixed a crash in 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).
• **Fixed 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.
• **C++ samplers no longer truncate the matrix-F degrees of freedom** to integers (non-integer delta from prior_sd).
• **Fixed 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.
• **Fixed 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.
• **Fixed qplot() deprecation warning**: Replaced deprecated qplot() with ggplot() + geom_density() in plot_prior().
• **Fixed 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.

Performance

• **Removed unused posterior arrays from the C++ samplers used by 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.

Maintenance

• Added rlang to Imports for proper use of .data pronoun in ggplot2 aesthetics.
• Added unit tests for multiple functions.
• **Refactored 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.
• **Fixed the 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.