                 Changes in version 2.2.0 (2026-09-24)                  

New features

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

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.

                 Changes in version 2.1.6 (2025-12-02)                  

Major changes to ordinal sampler

  - Stan-style latent centering for ordinal models: Both the Albert and
    Cowles ordinals samplers have been refactored to improve numerical
    stability and mixing in the presence of skewed or ceiling/floor
    ordinal items:
      - Thresholds are now initialized from the empirical category
        frequencies of each variable (i.e., cumulative proportions
        mapped to the probit scale), rather than arbitrary or equally
        spaced cut-points.
      - Latent variable draws (Z) are initialized at the expected value
        of the truncated normal (conditional mean), rather than uniform
        draws across the truncation interval. This ensures that ceiling
        or floor items start in the correct region of the latent space.
      - After each update/sweep of a latent column (Z_j), the column is
        recentered (mean ≈ 0) and the corresponding thresholds are
        shifted by the same offset, preserving the likelihood but
        aligning the latent origin. This mirrors the identification
        scheme used in Stan’s ordered-probit/ordered-logit models
        (latent mean fixed at 0, thresholds floating).
      - Probit‐scale bounds (±8 on the standard normal scale) have been
        introduced to cap semi-infinite truncation regions (e.g.,
        ((\tau_{K-1},∞))) to avoid numerical overflow and improve
        stability for extreme category distributions.

Bug fixes & improvements

  - Synchronized threshold matrices in the Cowles sampler: the
    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.
  - Improved sampler performance for ordinal variables with heavy tails
    or heavy ceiling/floor effects. This should lead to reduced
    shrinkage of partial-correlations toward zero when items are highly
    skewed.

Compatibility notes

  - The statistical model remains unchanged: you still get the same
    posterior for latent precision/correlation matrices. The changes are
    purely in parameterisation and initialization of the latent (Z) and
    threshold variables.

                 Changes in version 2.1.5 (2024-12-22)                  

Removed NPM library to avoid CRAN compiler errors

                 Changes in version 2.1.4 (2024-12-13)                  

Bug Fixes and Improvements

C++ search Function

1.  Improved Initialization
    
      - The initial adjacency matrix now takes start_adj (the maximum
        likelihood solution) as the starting point to avoid inefficient
        sampling.

2.  Adaptive Sampling
    
      - The sampling of zeros and nonzeros is now adapted to the newly
        accepted adjacency matrix (adj_s) rather than the static
        start_adj.

3.  Efficiency Enhancements
    
      - Skip updates on rejection since the graph remains unchanged.
      - Use find_ids(adj_mat) instead of find_ids(start_adj) to ensure
        edge modifications are correctly tracked after acceptance.

Minor Changes:

  - 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

                 Changes in version 2.1.3 (2024-07-05)                  

  - Replaced dprecated armadillo function conv_to<>::from with as_scalar
  - prior_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)

                 Changes in version 2.1.2 (2024-06-22)                  

  - The prior_sd (or rho_sd in var_estimate() ) is limited to ranges
    between 0 and sqrt(1/8). These values ensure that delta does not go
    below 1.
  - Critical: select() did not return partial correlations, but Fisher-z
    values in summary(). Fisher values are transformed back to
    correlation metric. This fixes #90, see changes.
  - Upgraded deprecated ggplot guides() argument
  - Resolved non positive definite initialization matrix in wishrnd() in
    copula models when NA's are present in observed variables (fixes
    #89). See changes here

                 Changes in version 2.1.1 (2024-02-23)                  

BFpack dependency error fixed.

                 Changes in version 2.0.1 (2020-07-23)                  

This version of BGGM included changes based on the JOSS reviews: see
here for the overview and here for specific issues.

                 Changes in version 2.0.0 (2020-05-31)                  

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.

Added features

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

Extensions

  - 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

                 Changes in version 1.0.0 (2020-02-06)                  

Initial CRAN release
