Fit a BayesR Mixture Model (MCMC or stochastic EM)

Description

BayesR (Erbe et al., 2012) places a four-component mixture prior on allele effects: zero, small, medium, and large effect classes with proportions pi_vec and class variances proportional to sigma2_ah. Use method = “mcmc” for full Bayesian posterior inference, or method = “em” for fast stochastic-EM point estimates (Gaussian only) (Wang et al., 2015).

Usage

run_bayesr(
  w,
  y,
  X = NULL,
  marker_type = c("auto", "snp", "multiallelic"),
  pi_vec = c(0.95, 0.02, 0.02, 0.01),
  sigma2_e_init,
  sigma2_ah = NULL,
  sigma2_vec = NULL,
  prior_params = NULL,
  mcmc_params = NULL,
  em_params = NULL,
  method = c("mcmc", "em"),
  response_type = c("gaussian", "binary"),
  fold_id = 0L,
  save_rds = FALSE,
  save_path = NULL,
  verbose = FALSE,
  map = NULL,
  windsize = NULL,
  windnum = NULL
)

Arguments

`w`	Numeric design matrix (`n x p`). Typically the `$W_ah` element returned by `construct_wah_matrix` or the `$W` element returned by `construct_snp_matrix`.
`y`	Phenotype vector (length `n`). For binary traits use 0/1 coding.
`X`	Optional fixed-effects design matrix (`n x q`). When supplied, the model is $y = X\alpha + W\beta + \mu + \epsilon$ with a flat prior on $\alpha$. Do not include a column of ones for the intercept; $\mu$ is sampled separately.
`marker_type`	One of `“auto”` (default), `“snp”`, or `“multiallelic”`. `“auto”` resolves to `“multiallelic”`. Selects the default prior for the BayesR sampler — see `sigma2_ah`, `sigma2_e_init`, `prior_params`. SNP path additionally rescales `variance_class` by `ncol(w) / sum(apply(w, 2, var))`.
`pi_vec`	Initial mixture proportions for the four classes (zero, small, medium, large). Must sum to 1. Default `c(0.95, 0.02, 0.02, 0.01)`.
`sigma2_e_init`	Initial residual variance. Default for `“snp”`: `var(y) / 2`. `“multiallelic”`: required (typically `var(y) * 0.5`; for binary use `1.0`).
`sigma2_ah`	Prior total genetic variance. Default for `“snp”`: `var(y) / 4`. `“multiallelic”`: required (typically `var(y) * 0.5`; for binary use `1.0`).
`sigma2_vec`	Optional explicit variance vector for EM; if `NULL` it is derived from `sigma2_ah`, `pi_vec`, and `prior_params$variance_class`.
`prior_params`	Optional named list overriding defaults. User- supplied fields always win. Recognised: `a0_e`, `a0_g`, `variance_class`. Defaults for `“snp”`: `a0_e = -1`, `a0_g = 2`, `variance_class = c(0, 0.0001, 0.001, 0.01)`. For `“multiallelic”`: `a0_e = 10`, `a0_g = 10`, `variance_class = c(0, 0.001, 0.01, 0.1)`. `b0_e` and `b0_g` are computed internally.
`mcmc_params`	Optional named list: `n_iter` (40000), `n_burn` (20000), `n_thin` (10), `seed` (123).
`em_params`	Optional named list: `max_iter` (500), `tol` (`1e-6`).
`method`	Either `“mcmc”` (default) or `“em”`.
`response_type`	`“gaussian”` (default) or `“binary”`. Binary requires `method = “mcmc”`.
`fold_id`	Integer label printed in progress messages; useful when running CV loops.
`save_rds`	If `TRUE` (default) the fit object is written to disk as an RDS file. Set `FALSE` (default) for CV loops.
`save_path`	Optional explicit RDS path. If `NULL` (default), defaults to `“results_bayesr.Rds”` in the current working directory.
`verbose`	If `TRUE`, the Rust engine streams per-iteration progress (start banner, Iter X/Y diagnostics, ESS / Geweke, completion) to stderr. Default `FALSE` keeps the Rust side silent. The brief R post-fit summary (model, runtime, $h^2$, ESS, Geweke) prints regardless.
`map`	Optional GWAS map. Supplying `map` (together with exactly one of `windsize` or `windnum`) computes per-allele PIP, per-block PIP, and per-window WPPA from the posterior samples and attaches them to the fit. For SNP markers, supply a data.frame with columns `SNP`, `CHROM`, `POS` (one row per W column). For multi-allelic (MH) markers, supply a data.frame with columns `block_id`, `chr`, `start_pos`, `end_pos` (one row per unique block_id in `attr(w, “block_id”)`); this matches the schema of `microhaplotype_coordinates.csv` from the `maspipeline` tool. Marker type is auto-detected from `attr(w, “block_id”)`. `map = NULL` (default) skips the GWAS branch and returns a fit object with no `pip`/`gwas` fields.
`windsize`	Window size in base pairs for WPPA. Mutually exclusive with `windnum`. Required when `map` is supplied unless `windnum` is.
`windnum`	Number of consecutive markers (or MH blocks) per window for WPPA. Mutually exclusive with `windsize`.

Details

Statistical model. The base linear model for BayesR is:

y = + W+ , (0,,_e^2 I)

where $y$ is the phenotype vector, $\mu$ is the intercept, $W$ is the marker design matrix (from construct_wah_matrix for multi-allelic or construct_snp_matrix for SNP markers), and $\beta$ is the allele-effect vector. When fixed effects are supplied via X:

$y = \mathbf{1}\mu + X\alpha + W\beta + \varepsilon$

BayesR mixture prior.

j ,,^2 ;; _0,(0,0) + 1,(0,,c_1^2) + 2,(0,,c_2^2) + 3,(0,,c_3^2)

where $(c_1,c_2,c_3)$ are the nonzero classes of variance_class and $\pi$ is updated each MCMC iteration via Dirichlet posterior. Residual: $\sigma_e^2 \sim \mathrm{InvGamma}(a_{0e},\,b_{0e})$. Genetic base scale: $\sigma_g^2 \sim \mathrm{InvGamma}(a_{0g},\,b_{0g})$.

Algorithm choice. MCMC uses marginalised Gibbs sampling and returns full posterior chains; recommended when posterior uncertainty, ESS, or Geweke diagnostics are needed. EM is much faster but yields no uncertainty quantification and does not support response_type = “binary”.

Auto-save. By default the returned fit is also saved to results_bayesr.Rds in getwd(). Set save_rds = FALSE (or pass an explicit save_path) when running cross-validation loops to avoid overwriting earlier folds.

Three usage scenarios. The same run_bayesr() + summary() + predict() pipeline supports:

Full-data fit: train on all data; predict(fit) returns in-sample metrics.
Train-test split: train on a subset; evaluate via predict(fit, W_test, y_test).
k-fold CV: the user loops over folds, calling run_bayesr() with save_rds = FALSE and fold_id = k, then predict() on each held-out fold.

Value

An object of class c(“masbayes_bayesr”, “masbayes”) — a list with the following key fields:

beta_hat, mu_hat, sigma2_e_hat: Posterior point estimates.
beta_samples, gamma_samples, pi_samples, sigma2_e_samples, sigma2_small_samples, sigma2_medium_samples, sigma2_large_samples, mu_samples: Posterior chains (single-row matrices for EM).
GEBV / pred_train: Training genomic estimated breeding values. For binary traits this is on the liability scale.
prob_train: Binary only: training-set predicted probabilities P(y = 1) = pnorm(pred_train) (probit inverse link from Albert-Chib augmentation).
h2, sigma2_g, sigma2_e: Heritability and variance components.
runtime: Elapsed seconds.
training_metrics: R2, RMSE, accuracy (or AUC for binary), bias. For binary, all metrics are on the observed (probability) scale — bias is the calibration slope (1.0 = perfectly calibrated) and RMSE^2 approximates the Brier score. AUC is rank-invariant so unaffected.
diagnostics: ESS and Geweke Z for key parameters (MCMC only).
variance_components: Posterior mean and 95% CI for each class plus the mixture proportions pi.
rds_path: Path of the saved RDS file, or NULL when save_rds = FALSE.
pip: (GWAS only, since v0.5.0) Per-allele posterior inclusion probability, length ncol(w). Computed as colMeans(gamma_samples != 0). NULL unless map was supplied.
pip_block: (GWAS only) Per-block PIP, length nrow(map). For SNP markers, identical to pip. For multi-allelic markers, a block is "active" in a posterior sample if at least one of its alleles has gamma != 0. NULL unless map was supplied.
gwas: (GWAS only) Data frame with columns Wind, Chr, N, Start, End, WPPA. Each row is one physical window; WPPA is the window posterior probability of association, i.e. the fraction of posterior samples in which at least one marker / block in the window has non-zero effect.
gwas_meta: (GWAS only) List capturing the call: windsize, windnum, midpoint_convention, marker_type, n_windows.

GWAS

Supplying a map together with exactly one of windsize or windnum turns run_bayesr() into a one-shot GWAS engine alongside genomic prediction. The Rust kernel is untouched; PIP and WPPA are derived from the existing gamma_samples matrix.

Map schema. For SNP markers, map is a data.frame with columns SNP (character), CHROM (integer), POS (integer), one row per W column. For multi-allelic (MH) markers, map has columns block_id (character), chr (integer), start_pos (integer), end_pos (integer), one row per unique block_id in attr(w, “block_id”). The MH schema matches the microhaplotype_coordinates.csv produced by the maspipeline tool, so production pipelines can pass it through without translation.

Windows. windsize (bp) closes a window when the next marker’s position exceeds the window start by windsize. windnum groups windnum consecutive markers / blocks per window. Windows never cross chromosome boundaries.

Constraints. GWAS requires method = “mcmc” (EM has no posterior samples). run_bayesa accepts the same arguments for API symmetry but errors out (BayesA has no zero-effect mixture mass, so PIP / WPPA are ill-defined). map = NULL (default) skips the GWAS branch entirely.

References

Erbe, M., Hayes, B. J., Matukumalli, L. K., Goswami, S., Bowman, P. J., Reich, C. M., Mason, B. A., & Goddard, M. E. (2012). Improving accuracy of genomic predictions within and between dairy cattle breeds with imputed high-density single nucleotide polymorphism panels. Journal of Dairy Science, 95(7), 4114-4129. doi:10.3168/jds.2011-5019

Wang, T., Chen, YP. P., & Goddard, M. E. (2015). A computationally efficient algorithm for genomic prediction using a Bayesian model. Genet Sel Evol, 47(34). doi:10.1186/s12711-014-0082-4

Examples

Code

library("masbayes")

d <- load_data("small")
mcmc <- list(n_iter = 1000L, n_burn = 500L, n_thin = 5L, seed = 123L)

# ---- (A) SNP path -----------------------------------------------------
# Pair `encoding = "zscore"` on the matrix builder with
# `marker_type = "snp"` on the fitter for the alternative biallelic-SNP
# convention.
snp_train <- construct_snp_matrix(d$snp[d$train_idx, ], encoding = "zscore")
W_train   <- snp_train$W
y_train   <- d$pheno$y_cont_qtl_snp[d$train_idx]
X_train   <- model.matrix(~ sex - 1, data = d$pheno[d$train_idx, ])

fit_snp <- run_bayesr(
  w             = W_train,
  y             = y_train,
  X             = X_train,
  marker_type   = "snp",
  sigma2_e_init = var(y_train) * 0.5,
  sigma2_ah     = var(y_train) * 0.5,
  mcmc_params   = mcmc,
  save_rds      = FALSE
)
summary(fit_snp)

# ---- (B) Microhaplotype path -----------------------------------------
# d$mh is consumable directly by construct_wah_matrix via attr("block_id").
# d$allele_freq is the frequency table required for the training call.
bid    <- attr(d$mh, "block_id")
hap_tr <- d$mh[d$train_idx, ]
W_mh   <- construct_wah_matrix(hap_tr, bid, d$allele_freq)$W_ah
y_mh   <- d$pheno$y_cont_qtl_mh[d$train_idx]

# ---- (C) GWAS (BayesR-only) ------------------------------------------
# Supply `map` together with `windsize` (bp) or `windnum` (markers per
# window) to compute per-allele PIP, per-block PIP, and per-window WPPA.
# Schema for SNP map: SNP / CHROM / POS (one row per W column).
# The map's row order does not need to be pre-sorted.
map_snp <- data.frame(
  SNP   = paste0("M", seq_len(ncol(W_train))),
  CHROM = rep(seq_len(4L), each = ncol(W_train) %/% 4L),
  POS   = rep(seq(1e6, by = 1e5,
                  length.out = ncol(W_train) %/% 4L), 4L)
)
fit_gwas <- run_bayesr(
  w             = W_train,
  y             = y_train,
  marker_type   = "snp",
  sigma2_e_init = var(y_train) * 0.5,
  sigma2_ah     = var(y_train) * 0.5,
  mcmc_params   = mcmc,
  save_rds      = FALSE,
  map           = map_snp,
  windsize      = 5e5
)
head(fit_gwas$gwas[order(-fit_gwas$gwas$WPPA), ], 5)

fit_mh <- run_bayesr(
  w             = W_mh,
  y             = y_mh,
  X             = X_train,
  sigma2_e_init = var(y_mh) * 0.5,
  sigma2_ah     = var(y_mh) * 0.5,
  mcmc_params   = mcmc,
  save_rds      = FALSE
)
summary(fit_mh)