Fit a BayesA Marker-Specific Variance Model (MCMC or stochastic EM)

Description

BayesA (Meuwissen et al., 2001) places a scaled inverse chi-squared prior on each allele’s effect variance, allowing heavier-tailed effect-size distributions than ridge regression. Use method = “mcmc” for full posterior inference or method = “em” for fast stochastic-EM point estimates (Gaussian only).

Usage

run_bayesa(
  w,
  y,
  X = NULL,
  marker_type = c("auto", "snp", "multiallelic"),
  nu = 4.5,
  sigma2_g,
  sigma2_e_init,
  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`). Use 0/1 for binary traits.
`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 BayesA sampler — see `nu`, `sigma2_g`, `sigma2_e_init`, `prior_params`.
`nu`	Degrees of freedom for the scaled inverse chi-squared prior on marker variances. Smaller values give heavier tails. Must be > 2. Default: `4` for `“snp”`, `4.5` for `“multiallelic”`.
`sigma2_g`	Prior total genetic variance. Default for `“snp”`: `var(y) / 4`. `“multiallelic”`: required (typically `var(y) * 0.5`; for binary use `1.0`).
`sigma2_e_init`	Initial residual variance. Default for `“snp”`: `var(y) / 2`. `“multiallelic”`: required (typically `var(y) * 0.5`; for binary use `1.0`).
`prior_params`	Optional named list. Only `a0_e` is used; `b0_e` is computed as `sigma2_e_init * (a0_e - 1)`, except when SNP-mode flat default applies (`a0_e = -1`, `b0_e = 0`). User-supplied `a0_e` always wins. Default `a0_e`: `-1` for `“snp”`, `10` for `“multiallelic”`.
`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 for progress messages.
`save_rds`	If `TRUE`, the fit object is saved to disk as an RDS file. Set `FALSE` (default) for CV loops.
`save_path`	Optional explicit RDS path. If `NULL` (default), defaults to `“results_bayesa.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 prints regardless.
`map`, `windsize`, `windnum`	Accepted for API symmetry with `run_bayesr` but not supported for BayesA. Passing any non-NULL `map` raises an error; BayesA has no zero-effect mixture mass, so PIP and WPPA are ill-defined for it. Use `run_bayesr` for GWAS-style posterior summaries.

Details

Statistical model. The base linear model for BayesA 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$

BayesA prior.

j {j}^2 ;; (0,,{j}^2), {_j}^2 ;; ^{-2}(,,S2)

Residual: $\sigma_e^2 \sim \mathrm{InvGamma}(a_{0e},\,b_{0e})$.

Algorithm choice. MCMC uses marginalised Gibbs sampling and returns full posterior chains. EM is much faster but provides no uncertainty quantification and does not support response_type = “binary”.

Auto-save. By default the returned fit is saved to results_bayesa.Rds in getwd(). Set save_rds = FALSE for CV loops.

Three usage scenarios. run_bayesa() + summary() + predict() support full-data fits, train-test splits, and k-fold CV uniformly. See run_bayesr for example loops.

Value

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

beta_hat, mu_hat, sigma2_e_hat, sigma2_j_hat: Posterior point estimates (intercept, allele effects, residual variance, per-marker variances).
beta_samples, sigma2_j_samples, sigma2_e_samples, mu_samples: Posterior chains (single-row matrices for EM).
GEBV / pred_train, h2, sigma2_g, sigma2_e: Training GEBVs (liability scale for binary), heritability, and total variance components.
prob_train: Binary only: training-set predicted probabilities P(y = 1) = pnorm(pred_train) (probit inverse link from Albert-Chib augmentation).
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 / Geweke Z (MCMC only).
variance_components: Per-marker variances binned into tertiles (small / medium / large) reporting mean, range, and marker count.
rds_path: Path of the saved RDS file, or NULL.

References

Meuwissen, T. H. E., Hayes, B. J., & Goddard, M. E. (2001). Prediction of total genetic value using genome-wide dense marker maps. Genetics, 157(4), 1819-1829. doi:10.1093/genetics/157.4.1819

Examples

Code

library("masbayes")

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

# ---- (A) SNP path -----------------------------------------------------
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_bayesa(
  w             = W_train,
  y             = y_train,
  X             = X_train,
  marker_type   = "snp",
  nu            = 4.5,
  sigma2_g      = var(y_train) * 0.5,
  sigma2_e_init = var(y_train) * 0.5,
  mcmc_params   = mcmc,
  save_rds      = FALSE
)
summary(fit_snp)

# ---- (B) Microhaplotype path -----------------------------------------
# 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]

fit_mh <- run_bayesa(
  w             = W_mh,
  y             = y_mh,
  X             = X_train,
  nu            = 4.5,
  sigma2_g      = var(y_mh) * 0.5,
  sigma2_e_init = var(y_mh) * 0.5,
  mcmc_params   = mcmc,
  save_rds      = FALSE
)
summary(fit_mh)