Preparing summaries, fitting Bayesian models and scoring

These interfaces connect existing native components. Preparation is explicit; estimators do not silently align alleles or change the marker universe. Current examples are shown below.

Prepare once, then choose the analysis

For already aligned input, no preparation call is necessary. Otherwise, gbase::checkStat() maps user columns to an ordinary qgg-style checked data frame, validates values, harmonizes against Glist alleles and applies explicit QC thresholds:

stat <- gbase::checkStat(Glist = Glist, stat = raw)

Here raw contains explicitly mapped marker, coordinate, allele and analysis columns. checkStat() preserves the ordinary data-frame contract with canonical names such as rsids, chr, pos, ea, nea, eaf, b, seb, stat, p, n and info. Exact and swapped alleles can be harmonized explicitly; ambiguous or incompatible records are rejected or removed according to the requested policy. It does not standardize effects or infer the generating model. Task packages validate the subset and scale they require.

Association evidence for marker sets

gbase::checkStat() can perform explicit mapping and Glist harmonisation; the association input remains an ordinary checked data frame. glma::vegas() calculates marker-set evidence against sparse LD attached by gbase::gprep(task = "sparseld"):

evidence <- glma::vegas(statistics, Glist, sets = marker_sets)
evidence$stat

Sets may represent genes, regions or annotations. This reuses native glma quadratic tests, with LD-aware null distributions and between-set score dependence. It is not enrichment or a claim of executable/RNG equivalence with historical VEGAS software. See the GWAS-to-pathway example for the connection to gsea() and its assumptions.

Construct and apply genomic scores

weights <- gscore::clump_threshold(statistics, reference_Glist,
  p_threshold = 5e-8, r2_threshold = 0.1, window_bp = 250000)
scores <- gscore::gscore(weights, target_Glist)
scores$scores

clump_threshold() is the only scoring-specific weight construction route. Its sparse-LD reference must cover the whole clumping rule: no variant-count window (max_distance_variants = 0), a base-pair window at least as wide as window_bp, and an LD construction r2 threshold no larger than r2_threshold. The bounded reference in the package tutorial is not suitable for its 250 kb C+T example; the scoring tutorial shows a separate reference construction. Other weights come from their statistical model owner, for example gbayes_fit$weights, including method = "sbayesrc". The portable weight frame declares marker coordinates, effect and other alleles, allele frequency, and one or more numeric score columns. gscore() aligns exact or swapped alleles to every active Glist BED component, preserves selected sample order, and records omitted weighted markers. It supports standardized or raw mean-imputed dosage coding and applies multiple score columns in one call.

Bayesian posterior means can be applied genome-wide or by supplied marker sets:

fit <- gbayes::gbayes(
  checked, reference_Glist, method = "bayesr",
  control = list(burnin = 1000, sampling_sweeps = 5000,
                 seeds = c(11, 29))
)

scores <- gscore::gscore(
  fit$weights, target_Glist,
  sets = marker_sets,
  include_total = TRUE
)

sets is a uniquely named list of marker names. The target Glist resolves each name to its BED component and physical column, so callers do not split sets by chromosome or calculate file offsets. Sets may overlap and span BED files. The native scorer decodes each required marker once per file and accumulates all set scores from sparse membership, without constructing a dense marker-by-set coefficient matrix. scores$coverage, scores$set_membership and scores$missing_set_markers record the resolved contract. By default only set scores are returned; include_total=TRUE also retains the original score columns. Optional named set_weights are applied after score standardization.

gscore does not fit ridge, L1 or mixture models and does not select models on validation summaries. Those are statistical-learning tasks, not genomic scoring.

Additional Bayesian inference engines

checked <- gbase::checkStat(Glist=Glist, stat=reported)
vb <- gbayes::gbayes(checked, Glist, method="variational",
                     prior=list(p=.01, h2=.4))
mixture <- gbayes::gbayes(checked, Glist, method="mixture_variational",
                  prior=list(variances=c(0,.01), weights=c(.9,.1)))
gibbs <- gbayes::gbayes(checked, Glist, method="point_normal",
                prior=list(p=.01, h2=.4),
                control=list(seed=17, burnin=500, sampling_sweeps=1000))
summary(vb)

These routes accept one trait and a Glist containing prepared sparse LD, and return the common gbayes_fit contract with method-specific estimates and diagnostics.

  • Point-normal Gibbs: p and h2 determine the initial prior. Controls include burnin=500, sampling_sweeps=200, seed=1, sparse=FALSE and auto_update=FALSE. Native parameter paths remain visible. Marker variances are not supplied by this engine and remain unavailable.
  • Point-normal / mixture variational: fixed priors, maximum_iterations=200 and tolerance=1e-8. Optional windows is a two-column matrix of one-based inclusive ranges. Windows may overlap; a single ordered window pass does not establish whole-panel convergence. Objective/drift/window diagnostics remain available; posterior variances are variational approximations.
  • Continuous shrinkage: prior=list(global_scale=1); controls burnin=500, sampling_sweeps=500, thinning=5, seed=1, learn_global_scale=FALSE, maximum_block_dimension=2048 and required block_offsets=c(0,...,M). This is an explicitly zero-based boundary vector, not marker IDs. Native checks reject retained cross-block LD. LD must be materialized. Local/global/residual-scale summaries are returned; this continuous prior has no inclusion PIP, and marker variances remain unavailable.

These bindings use one chain/fit (threads=1). They do not inherit BayesC/BayesR multitrait, annotation, parameter-learning or posterior-summary controls. They do not manufacture ESS, R-hat or uncertainty absent from the native result. Independent reruns and appropriate statistical qualification remain necessary.

Reuse GNOVA annotations for grouped covariance models

library(gcorr)
grouped <- gcorr_grouped_covariance(source_fit, continuous = TRUE,
  units = setNames(rep("standardized genetic covariance", 3),
                   c("A", "B", "C")))
sem <- gcorr_grouped_sem(grouped, sem_model)
ggm <- gcorr_grouped_ggm(grouped, edges)

gcorr accepts a coordinated continuous-trait GNOVA or SumHer fit that retains full cross-group jackknife covariance. It returns joint SEM or fixed-graph GGM estimates and exposes explicit cross-group test functions. See the covariance guide and the deterministic example. Regional LAVA partial-correlation and regression analyses remain separate workflows.