Fit BayesC and BayesR to the simulated individual-level panel and inspect posterior parameter estimates. We use all individuals; prediction and the historical eight-fit comparison are in the advanced example. Use gbayes 0.1.5 or later and a compatible gbase installation. This page shows a runnable recipe; rendering does not fit the models.
Set the working directory to an existing folder of your choice. Download the same panel, corrected marker map and phenotype used in the linear and mixed models example.
library(gbase)
library(gbayes)
setwd("C:/Users/Projects/Examples")
url <- "https://github.com/psoerensen/qgdata/raw/main/simulated_human_data/"
base <- "https://psoerensen.github.io/gact/Document/glma-simulated-data/"
download.file(paste0(url, "human.bed"), "human.bed", mode="wb")
download.file(paste0(base, "human.bim"), "human.bim")
download.file(paste0(url, "human.fam"), "human.fam")
Glist <- gprep(study="Simulation", bedfiles="human.bed",
bimfiles="human.bim", famfiles="human.fam")
rsids <- gfilter(Glist=Glist, excludeMAF=0.05, excludeMISS=0.05,
excludeCGAT=TRUE, excludeINDEL=TRUE, excludeDUPS=TRUE,
excludeHWE=1e-12, excludeMHC=FALSE)
Glist <- gprep(study="Simulation", bedfiles="human.bed", rsids=rsids)
phenotypes <- read.csv(paste0(base, "example-phenotypes.csv"),
colClasses=c("character", "numeric"))
y <- setNames(phenotypes$phenotype, phenotypes$id)
data <- data.frame(id=Glist$ids, phenotype=as.numeric(y[Glist$ids]))fit_c <- gbayes(data=data, Glist=Glist, trait="phenotype", method="bayesc")
fit_cThe defaults learn effect variance, residual variance and mixture weights. They use two reproducible chains, 1000 burn-in and 2000 sampling sweeps, thinning by five and full marker sweeps. Starting variances use the observed phenotype variance; they do not use the known number of causal markers. The phenotype is centred internally without rescaling, so estimates retain its original units. These are starting settings, not a convergence guarantee.
fit_c$posterior
head(fit_c$estimates)| Parameter | Meaning |
|---|---|
VB (effect_variance) |
Slab variance scale for standardized marker effects |
VE (residual_variance) |
Residual variance |
VG (genetic_variance) |
Realized genetic variance in this panel |
VLE (vle) |
Sum of individual-marker variance contributions |
VLD (vld) |
Sum of cross-marker contributions, VG − VLE |
h² (h2) |
Draw-wise VG / (VG + VE) |
The table reports posterior means, SDs, credible intervals, R-hat, bulk ESS and Monte Carlo error. VLD can be negative. Posterior mean VG also differs from the variance of scores calculated using posterior mean marker effects.
fit_r <- gbayes(data=data, Glist=Glist, trait="phenotype", method="bayesr")
fit_r
fit_r$posteriorBayesR uses several nonzero variance components. Its initial mixture
weights and its Dirichlet learning prior are separate settings. The
fitted $prior and $control show the resolved
defaults; explicit settings remain available. See the defaults
guide for the initialization formulas and explicit-prior
compatibility.
fit_c$posterior[, c("parameter", "rhat", "ess_bulk", "mcse_mean")]
fit_r$posterior[, c("parameter", "rhat", "ess_bulk", "mcse_mean")]Printing a fit flags estimated quantities with unavailable diagnostics, R-hat above 1.01 or bulk ESS below 100. Inspect traces and continue sampling when needed; passing this screen does not prove convergence. Check variance and mixture quantities separately. The advanced investigation shows why stable predictions can coexist with poor parameter mixing.
fit_scheduled <- gbayes(data=data, Glist=Glist, trait="phenotype",
method="bayesc",
control=list(marker_schedule=list(enabled=TRUE, burnin_only=TRUE)))
fit_scheduled$posteriorPosterior sampling still visits every marker. Compare diagnostics before interpreting differences between fits. The simulated independent-marker panel illustrates the variance decomposition; it does not validate realistic population LD.