Bayesian Alphabet TWAS Weights From Summary Statistics — bayesAlphabetRssWeights • pecotmr

Fits a Bayesian alphabet prior to marginal summary statistics and an LD reference and returns the posterior mean effects. Each method runs on several engines, chosen by backend:

"qgg": qgg::gbayes on the LD matrix, for bayesA, bayesC, bayesL, bayesN and bayesR.
"hibayes": hibayes::sbrm on the LD matrix, for every method.
"pecotmr": pecotmr's own sampler, on a genotype reference's factor of the LD without forming the matrix. Its priors follow qgg's for bayesA, bayesC, bayesN and bayesR and hibayes' for bayesB; bayesL takes Park and Casella's form (hibayes') or Legarra et al.'s (qgg's individual-level one), through BayesAlphabetRssOptions. By default it reads the statistics through SDPR's likelihood (Zhou and Zhao 2021), which drops the part of the marginal effects the reference LD cannot explain; BayesAlphabetRssOptions(ldProjection = NULL) turns that off (see ldProjection of the internal bayesAlphabetRss()).

NULL picks "pecotmr" for a genotype reference, and for an LD matrix qgg (bayesA, bayesC, bayesN, bayesR) or hibayes (bayesB, bayesL). qgg and hibayes run on a genotype reference through the LD matrix formed from it; pecotmr's sampler does not run on an LD matrix.

Usage

bayesAlphabetRssWeights(
  stat,
  ldReference,
  method,
  backend = NULL,
  pi = NULL,
  nIter = 5000,
  nBurnin = 1000,
  thin = 5,
  shrinkCorr = 0.95,
  fitRetention = c("none", "slim", "full"),
  methodArgs = NULL
)

Arguments

stat: A list of summary statistics: n (per-variant sample sizes) and the marginal effects as cor (standardized), z, b with seb, or b alone, taken as already standardized. The engines receive the standardized effects z / sqrt(n), with n the median of n.
ldReference: The LD reference the statistics' variants are matched to by stat$variantNames: a genotype panel (see readGenotypes), a GenotypeMixture or an LdMatrix.
method: The prior: "bayesA", "bayesB", "bayesC", "bayesL", "bayesN" or "bayesR".
backend: The engine, one the method has (see Description); NULL picks the default for the LD reference.
pi: For bayesB and bayesC, the prior proportion of non-null effects, where each engine's estimate of it starts; NULL gives 0.2 for bayesB and 0.1 for bayesC. Refused for the other methods.
nIter: Integer. Total number of MCMC iterations, burn-in included, on every engine. Default 5000.
nBurnin: Integer. Number of burn-in iterations discarded; less than nIter. Default 1000.
thin: Integer. Thinning interval for retained MCMC samples. Default 5.
shrinkCorr: The LD shrinkage $c$: every engine reads $cR + (1 - c)I$, qgg and hibayes as the matrix they are given, pecotmr's sampler through the factor. Default 0.95, the LDpred2 vignette's. A reference with fewer samples than the window has variants lacks directions the statistics have, and without shrinkage the samplers' effects can grow without bound along them (pecotmr's sampler also drops those directions by default, through SDPR's likelihood; see bayesAlphabetRssWeights). A small or mismatched reference does better with less, down to 0.4 as the LDpred2 vignette advises; the study's own genotypes as the reference need none, up to 1.
fitRetention: How much of the fit is kept as the "fit" attribute of the returned weights: "none" (default) attaches nothing, "slim" and "full" attach the engine's estimates (betaEst, and where the engine gives them the inclusion probabilities pip, the non-null proportion pi and the residual variance), backend and the shrinkCorr used. No engine keeps trimmable intermediates, so the two retaining levels behave alike.
methodArgs: Options for the backend's engine, built with its constructor: QggOptions for qgg, HibayesRssOptions for hibayes, BayesAlphabetRssOptions for pecotmr's sampler. NULL keeps the engine's defaults.

Value

A numeric vector of the posterior mean effects on the standardized scale, one per variant. Unless fitRetention = "none", the engine's estimates are attached as attribute "fit".

Examples

data(eqtlRegionExample)
X <- eqtlRegionExample$X[, 1:30]
y <- eqtlRegionExample$yRes
ss <- lapply(
  seq_len(ncol(X)), function(j) coef(summary(lm(y ~ X[, j])))[2, 1:2])
stat <- list(
  b = vapply(ss, `[`, numeric(1), 1L),
  seb = vapply(ss, `[`, numeric(1), 2L),
  n = rep(nrow(X), ncol(X)),
  variantNames = colnames(X)
)
w <- bayesAlphabetRssWeights(
  stat,
  LdMatrix(cor(X)),
  "bayesN",
  backend = "hibayes",
  nIter = 1000,
  nBurnin = 200
)