Construct a SNP Design Matrix

Description

Build a SNP design matrix from an allele-dosage matrix $X$ (entries 0, 1, 2). Two encodings are supported: VanRaden centering (default, $W_{ij} = X_{ij} - 2 p_j$) and per-column z-score standardisation ($W_{ij} = (X_{ij} - 2 p_j) / s_j$). The output is suitable for direct use with run_bayesr and run_bayesa.

Usage

construct_snp_matrix(
  X,
  encoding = c("vanRaden", "zscore"),
  ref_freq = NULL,
  ref_sd = NULL
)

Arguments

`X`	Numeric or integer matrix of allele dosages (typically `0`, `1`, or `2`) with dimensions `n x p_snp`. Coerced to `double` internally.
`encoding`	Either `“vanRaden”` (default) or `“zscore”`. Controls whether columns are additionally divided by their standard deviation after centering.
`ref_freq`	Optional numeric vector of length `ncol(X)` giving the reference allele frequencies $p_j$ from the training set. If `NULL` (default) frequencies are computed from `X` itself (training-set use).
`ref_sd`	Optional numeric vector of length `ncol(X)` giving the training-set column standard deviations. Required only when `encoding = “zscore”` and `X` is a test set. If `NULL` for a training call, the sd is computed from `X`. Ignored when `encoding = “vanRaden”`.

Details

Choosing an encoding.

“vanRaden” (default) — centers each column on the training-set allele frequency. Column variance is proportional to $2 p_j (1 - p_j)$, so rare variants contribute proportionally less to the implied genomic relationship matrix. This is the classical genomic-prediction convention and is consistent with tcrossprod(W) / k_grm used by GBLUP backends.
“zscore” — additionally divides each column by its standard deviation, so every marker contributes equal variance regardless of MAF. Aligned with alternative biallelic-SNP parameterisations that assume marker-uniform variance contributions.

Both encodings work with either marker_type = “multiallelic” or “snp” in the Bayesian fitters, but “zscore” is intended to pair with marker_type = “snp” for a fully alternative SNP convention.

Training / test workflow. Always centre (and scale, for zscore) the test set with training statistics, never with statistics recomputed from the test set itself:

Training: call construct_snp_matrix(X_train, encoding = …). The returned $freq and, for zscore, $sd hold the training statistics.
Test: call construct_snp_matrix(X_test, encoding = …, ref_freq = train$freq, ref_sd = train$sd) so test columns align with training columns under the same baseline.

For phased multi-allelic haplotype data, use construct_wah_matrix instead.

Value

A list with elements:

W: Numeric matrix n x p_snp. For vanRaden: $W_{ij} = X_{ij} - 2 p_j$. For zscore: $W_{ij} = (X_{ij} - 2 p_j) / s_j$.
freq: Numeric vector of length p_snp with the allele frequencies used for centering.
sd: Numeric vector of length p_snp with the column standard deviations used for scaling. Present only when encoding = “zscore”.
n, p: Number of individuals and SNPs.

Examples

Code

library("masbayes")

d <- load_data("small")
X_train <- d$snp[d$train_idx, ]
X_test  <- d$snp[d$test_idx, ]

# vanRaden (default)
train_v <- construct_snp_matrix(X_train)
test_v  <- construct_snp_matrix(X_test, ref_freq = train_v$freq)

# zscore (alternative SNP convention)
train_z <- construct_snp_matrix(X_train, encoding = "zscore")
test_z  <- construct_snp_matrix(X_test, encoding = "zscore",
                                ref_freq = train_z$freq,
                                ref_sd   = train_z$sd)

stopifnot(ncol(train_z$W) == ncol(test_z$W))