Bayesian genomic regression
The main interface is gbayes::gbayes():
fit <- gbayes::gbayes(
stat = NULL,
Glist = Glist,
data = NULL,
method = "bayesr",
prior = list(),
control = list(),
sampling = NULL,
annotation = NULL,
posterior = list(),
continuation = NULL
)stat + Glist selects a summary-statistics likelihood. data + Glist selects an individual-level BED likelihood. A named list of checked summary frames or several named phenotype columns selects the multivariate route. Unsupported method/data combinations fail explicitly.
Fixed-variance single-effect summaries
gbayes::gbayes_single_effect(stat) fits a deterministic single-effect model from an ordinary data frame with unique rsids and either b/seb or score/information. It accepts fixed positive residual_variance and effect_variance, optional nonnegative marker_weights, and output labels trait and response_units. No Glist or prepared LD resource is required by this focused interface.
The per-marker result contains posterior assignment probabilities and moments. The model assumes exactly one active marker; assignment probabilities sum to one even when the evidence is weak. They are not unconditional association probabilities. This interface does not expose native single-effect variance learning or the wider BayesC/BayesR control set. For regional multi-effect fitting and posterior interpretation, use the mapping guide.
Prepare genotypes and summaries
Single-trait summary BayesC/BayesR also accepts familiar iteration and prior arguments:
fit <- gbayes::gbayes(Glist=Glist, stat=stat, method="bayesC",
algorithm="mcmc", nit=100, nburn=0, seed=10)
head(fit$stat) # bm: native posterior mean; dm: PIP
head(fit$post) # ve, vg, vb, pi, pip on the fitted likelihood scale
head(fit$conv) # native R-hat, ESS, MCSE and diagnostic statusesFamiliar calls default to one chain, 100 sweeps, zero burn-in, h2=.5, BayesC active probability .001, and learned slab/residual variances and mixture weights. nit/nburn/nthin/nrun/seed, pi/gamma/vb/ve, and updateB/updateE/updatePi adapt to native settings. Explicit conflicting aliases and nested settings are rejected; existing advanced calls retain their defaults. Familiar adaptation is limited to single-trait summaries.
Summary bm remains in standardized native genotype/phenotype units. The standalone reference-LD residual policy defaults to adjustment .9; it differs from qgg’s marker-specific adjustment. No independent genetic-variance draw is recreated. post$vg instead averages the joint-effect quadratic b' C b / n from retained draws, reusing the fitted residual. It is unavailable with variable sample sizes or a distinct explicit posterior reference. pi is total non-null mixture probability and pip is summed marker PIP (expected active count); pi0, pi1, … retain component weights. conv exposes rhat_*, ess_*, mcse_* and status_* columns, rather than qgg Geweke z scores. Original native posterior and diagnostic fields remain available. Single-chain output does not establish convergence or posterior equivalence. Use gmap::getCredibleSets(Glist, fit, sets) to interpret saved marginal PIPs with full BED LD without refitting.
Scalar summary BayesC/BayesR support optional LD swaps on sparse-LD and block-eigen routes, with the same controls as regional gmap():
fit <- gbayes::gbayes(Glist=Glist, stat=stat, method="bayesC",
algorithm="mcmc", nit=1000, nburn=100, nrun=4, seed=10,
control=list(ld_swap=list(probability=.05, minimum_r2=.8,
maximum_friends=50L, moves=1L)))
fit$ld_swapNative Metropolis-Hastings proposals relocate an active effect to an inactive LD neighbour. probability selects sweeps for up to moves proposals; minimum_r2 and maximum_friends bound candidate neighbours. Neighbours use the fitted cross-products, including retained eigen truncation. fit$ld_swap reports pooled attempted/accepted counts and the acceptance rate; fit$chains retains per-chain counts. Counts include burn-in and no attempts gives an NA rate. They measure proposal activity, not convergence. Individual BED fits do not support these swaps. For sparse LD, neighbour preparation scans the existing stored correlations once and shares the list across chains. It does not recompute LD from genotypes or rebuild the list at every sweep.
Probability defaults to zero; zero moves also disables swaps and neighbour preparation. Finding neighbours adds preparation cost, especially with eigen operators. Existing manuscript timings used swaps disabled and remain unchanged; improved mixing or elapsed time has not been established by those timings.
BayesC and BayesR also accept algorithm="mcmc-eigen". Declare exhaustive contiguous LD blocks using gbase::define_ld_blocks(Glist, blocks); block edges must have no retained cross-block LD. The verified decomposition is prepared once per block and reused by every chain. This route uses the retained block-diagonal summary likelihood, defaults to positive eigenvalue mass .995, and learns a global residual variance in the retained response/rank space:
eigen_fit <- gbayes::gbayes(stat, Glist, method="bayesR",
algorithm="mcmc-eigen", nit=1000, nburn=100, nrun=4, seed=10,
control=list(eigen_control=list(retained_positive_mass=.995)))
eigen_fit$diagnostics$block_eigenUse updateE=FALSE for fixed residual variance. Reference-LD adjustment and caller-supplied full-space normalization cannot be combined with this route. SBayesRC remains a separate annotated model with block-specific residual scales. The eigen option is not a promise of lower elapsed time: preparation and retained rank affect the cost. The existing manuscript timings describe the earlier dense/sparse MCMC version and have not been replaced by these interface changes.
gbase::gprep() owns the qgg-style Glist contract. Sparse LD is added to the same object rather than creating another public descriptor.
Glist <- gbase::gprep(
bedfiles = c("chr1.bed", "chr2.bed")
)
Glist <- gbase::gprep(
Glist = Glist,
task = "sparseld",
out_prefix = "reference/ld"
)
stat <- gbase::checkStat(Glist = Glist, stat = reported_statistics)checkStat() returns an ordinary data frame with qgg-style names. It performs marker, coordinate and allele alignment but does not introduce a summary class, infer sample sizes or manufacture study frequencies. For the Gaussian summary likelihood, gbayes() requires finite b, seb and positive n; the package performs the documented standardized-effect conversion internally.
Scalar BayesC and BayesR
fit <- gbayes::gbayes(
stat,
Glist,
method = "bayesr",
prior = list(
residual_variance = 1,
effect_variance = 0.05,
weights = c(0.8, 0.15, 0.05),
variance_multipliers = c(0, 0.1, 1)
),
control = list(
burnin = 500,
sampling_sweeps = 1000,
seeds = c(11, 29)
),
trait = "height"
)For individual-level fitting, phenotype rows are matched by ID and may be in a different order from the selected Glist samples.
fit <- gbayes::gbayes(
data = phenotypes,
Glist = Glist,
trait = "height",
method = "bayesc",
prior = list(
residual_variance = 1,
effect_variance = 0.05,
inclusion_probability = 0.01
)
)Both routes support fixed or learned effect variance and mixture weights, explicit residual policies, compact posterior summaries, annotations and exact serialized continuation. Individual fits use native BED coding and return portable allele-labelled fit$weights.
Scheduled scalar BED fits and variance decomposition
With gbayes 0.1.2 or later, scalar individual BayesC and BayesR accept control$marker_schedule=list(enabled=TRUE). Defaults are full_sweep_every=10, base_interval=50, maximum_interval=200, candidate_threshold=.001, candidate_lifetime=20, and burnin_only=FALSE. A zero full-sweep interval requests full sweeps. Active and recent candidate markers are revisited; low-probability null markers can be skipped between full sweeps. Scheduling is an explicit approximation to ordinary full-sweep sampling. It is opt-in, scalar-BED-only in R, and its queue/counters are serialized for continuation.
For a posterior-analysis example, prefer scheduling only during burn-in:
control <- list(burnin = 1000, sampling_sweeps = 2000,
seeds = c(17, 29), retain_diagnostics = TRUE,
marker_schedule = list(enabled = TRUE, burnin_only = TRUE))These counts are starting settings, not a convergence guarantee. With burnin_only=TRUE, every sampling sweep is full, including continuation batches with zero burn-in. With scheduling during sampling, learned global mixture weights update on periodic full sweeps and, in corrected scalar native builds, whenever the scheduled visit includes every marker. Visits that skip markers hold weights between full sweeps; an interval of ten can therefore give ten times fewer weight updates. Earlier native builds also held weights on nonperiodic visits that included every marker, needlessly slowing learning in the all-hot case. A low candidate threshold can keep all markers hot, providing no marker-visit saving. The all-hot correction changes scheduled draw sequences; continue a study with a consistent native build.
Check posterior quantities separately: similar held-out accuracy and highly correlated posterior-mean predictions can coexist with poor variance or mixture-weight mixing. Inspect parameter-level R-hat, ESS, Monte Carlo error and chain traces; a single worst-parameter statistic does not describe every output. Starting weights or inclusion_probability are distinct from the Dirichlet/Beta weight_prior used when learning is enabled. Continuing a saved fit keeps its model and transition settings; changing a prior or scheduler requires a fresh fit. An extra continuation batch can be discarded as further warm-up before summarizing a subsequent batch. Full sweeps establish the baseline before studying approximate scheduling during sampling.
fit$individual_bed$variance exposes native final and post-burnin means:
[ V_g=|Xb|^2/n, V_{}=j(x_jx_j)b_j2/n, V{}=V_g-V_{}. ]
The LD term can be negative. Empirical diagonal variances need not be one under the native allele-frequency/HWE coding. These quantities differ from the slab variance parameter Vb; BayesR component variances are its multiplier-scaled values. Residual variance Ve is a separate model parameter. Compute h2 = Vg/(Vg+Ve) at each posterior draw when reporting model-scale heritability. The variance of posterior mean predictions differs from mean posterior Vg.
Each fit$chains[[k]]$individual_bed contains vb_trace, ve_trace, vg_trace, vle_trace, vld_trace, weight_trace, and log-CPO. BED traces include burn-in and every sweep; variance means average all post-burnin sweeps independently of effect thinning. Top-level traces average chains at corresponding sweeps; use per-chain traces for posterior uncertainty and diagnostics. Top-level final_* values average endpoints, not posterior draws. Per-chain scheduling counters are scheduled_marker_updates and scheduled_full_sweeps. Log-CPO uses retained training draws and is distinct from held-out prediction accuracy.
The gact individual-level example shows sparse-trait BayesC/BayesR with training-only preparation and held-out scoring. Its advanced study preserves the original short-run sparse/polygenic comparison and records quantity-specific prediction, variance and mixture diagnostics.
Annotation-informed block-eigen regression
method = "sbayesrc" fits a single checked summary trait with an annotation-informed BayesR prior. The Glist must contain a verified independent block partition. Each LD block is loaded and decomposed once per fit, and the prepared immutable factors are shared by all chains.
Glist <- gbase::define_ld_blocks(
Glist,
data.frame(
block_id = block_table$block_id,
first_marker = block_table$first_marker,
last_marker = block_table$last_marker
)
)
component_design <- as.matrix(annotation_table[fitted_markers, , drop = FALSE])
fit <- gbayes::gbayes(
stat,
Glist,
method = "sbayesrc",
annotation = list(component = list(design = component_design)),
control = list(
eigen_control = list(retained_positive_mass = 0.995),
block_eigen_residual = list(mode = "all_mix")
)
)The block residual update is defined on retained response coordinates: each block uses the squared retained response norm and retained rank. This is an explicit reduced-space likelihood and is not an assertion of identical normalization to GCTB or another SBayesRC implementation. Set retained_positive_mass = 1 to retain every positive eigenvalue. On a full-rank block-diagonal reference with fixed block residuals, this route reduces to the ordinary bayesr sufficient-statistic fit.
fit$diagnostics$block_eigen records block loads, eigensystem decompositions, retained ranks, discarded positive mass and peak dense-block memory. The factors remain native working state; the returned fit is serializable and continuation restores sampler state without serializing dense eigenvectors.
Additional summary engines
The same dispatcher provides focused single-trait summary methods:
point_normal <- gbayes::gbayes(stat, Glist, method = "point_normal")
shrinkage <- gbayes::gbayes(stat, Glist, method = "continuous_shrinkage")
variational <- gbayes::gbayes(stat, Glist, method = "variational")
mixture <- gbayes::gbayes(stat, Glist, method = "mixture_variational")These retain method-specific controls and outputs. They do not silently inherit all BayesC/BayesR sampling, learning or annotation features.
Multivariate and study-specific fits
Named phenotype columns fit a joint individual-level model:
joint <- gbayes::gbayes(
data = phenotypes,
Glist = Glist,
trait = c("height", "weight"),
method = "bayesr",
prior = list(
effect_covariance = effect_covariance,
residual_covariance = residual_covariance
)
)For summaries, stat may be a named trait list with one shared Glist sparse-LD reference, or a named study list paired with a named list of study-specific Glist references. Study identity remains distinct from trait identity. Local LD operators are preserved; no synthetic union LD matrix is created. Independent and shared-LD sampling assumptions are declared through sampling.
Annotations and continuation
BayesC supports inclusion and variance annotations. BayesR supports component and variance annotations. BayesS scaling uses the same checked variance-model boundary. Multivariate fits additionally support pattern annotations.
Pass a previous fit, or its $continuation member, back through continuation. The model-defining inputs, marker order, prior, annotations and transition controls must remain identical. Continuation restores native chain states and performs zero additional burn-in.
Posterior genetic quantities
An explicitly aligned sparse-LD Glist can define the reference population for posterior genetic covariance and correlation. Heritability additionally requires an explicit phenotype variance; no denominator is inferred.
groups <- setNames(marker_groups, fitted_marker_ids)
fit <- gbayes::gbayes(
stat,
Glist,
method = "bayesr",
posterior = list(
reference = Glist,
partition = groups,
phenotype_variance = c(height = 1)
)
)fit$quantities separates genetic covariance/correlation, heritability, reporting partitions, predictions and per-chain quantities. Partition groups are reporting objects and are independent of covariance-group priors.
Predict selected individuals with uncertainty
Posterior score uncertainty uses an existing target Glist. Its active markers may be a strict subset of fitted markers and may span any number of BED components. Components and markers retain Glist order; requested individuals retain the supplied ID order.
target <- gbase::gprep(
bedfiles = c("target-chr1.bed", "target-chr2.bed"),
ids = selected_ids,
rsids = selected_markers
)
fit <- gbayes::gbayes(
stat,
Glist,
method = "bayesr",
posterior = list(
prediction = list(
Glist = target,
ids = selected_ids,
target_set = "validation",
retain_draws = FALSE
)
)
)The package validates exact fitted-marker availability, coordinates and counted alleles. Component projections use the public gscore::summary interface and are summed into one score. Coverage records fitted, used and omitted markers and the number of active BED components. Reported score covariance is posterior effect uncertainty, not phenotype predictive variance.
Scalar individual-level BED fits also include genetic_variance (VG), vle, vld and h2 directly in fit$posterior. These rows use the original per-chain draws from every post-burn-in sweep, independently of effect thinning. They report posterior mean, SD, credible intervals, R-hat, ESS and Monte Carlo error using the native summary routines. h2 is summarized from each draw’s VG / (VG + VE). Existing effect_variance (VB) and residual_variance (VE) rows retain their effect-thinned parameter draws. Continuation summaries cover the new batch.
Apply portable posterior weights
When genotype coding is identified, fit$weights is directly consumable by gscore(). This includes block-eigen SBayesRC fits; scoring uses the portable posterior mean weights and does not need the fitted eigensystems:
prediction <- gscore::gscore(fit$weights, target)For pathway or other marker-set scores, pass a uniquely named list of marker IDs. The target Glist resolves their BED locations and the native scorer handles overlapping or cross-file sets in one pass:
prediction <- gscore::gscore(
fit$weights, target,
sets = marker_sets,
include_total = TRUE
)Weight learning remains with gbayes; score application remains with gscore. Raw association effects are not automatically Bayesian weights, and conventional clumping and thresholding remains gscore::clump_threshold().
Regional interpretation
gmap::gmap_bed() can consume a saved individual-BED gbayes_fit together with sparse LD prepared from the same Glist samples and physical marker selection:
mapped <- gmap::gmap_bed(fit, Glist, regions)This interprets global marginal PIPs by region without refitting. Marginal-PIP groups are not component credible sets.
gmap_bed() is specific to individual-BED fits because it verifies the fitted sample selection, BED files and in-sample LD identity. A summary SBayesRC fit remains a whole-reference gbayes_fit; regional BayesC/BayesR or multi-effect models are fitted through gmap() with their own explicit regional likelihood.
Qualification status
The package passed its complete source check with no errors or warnings, focused deterministic tests, an exact Windows DLL-closure and SHA-256 audit, and clean binary loading. The two source-check NOTEs concern installed DLL size and the conservative compiled-code scan of the bundled native DLL closure.
Simple individual-level defaults
Scalar individual BED fits can omit prior: gbayes initializes half the observed phenotype variance as residual variance and half as nominal genetic variance. BayesC starts with inclusion probability .001; BayesR starts with weights (.95,.02,.02,.01) and multipliers (0,.01,.1,1). The marker-scale variance is VG/(m times mixture variance mass), using standardized genotype coding. This is a starting heuristic, not an empirical estimate of heritability or causal count. Variance learning uses df=4 and scale=initial variance/2; mixture learning uses unit Beta/Dirichlet concentrations. These learning priors are separate from initial mixture probabilities. Automatic runs use 1000 burn-in and 2000 sampling sweeps, thinning 5, seeds 1 and 2, one thread, retained diagnostics and full sweeps. Control overrides remain available. Explicit nonempty priors retain the existing opt-in learning controls and require proper learning priors when enabled. Scalar BED responses are centred internally without rescaling; the default likelihood dimension is n-1. The phenotype mean is recorded in provenance. An omitted prior during continuation reuses the saved prior and transition controls. Inspect $prior, $control and $posterior; default sweep counts and the R-hat/ESS screening message do not guarantee convergence.