LDpred2-auto Weights from Summary Statistics

LDpred2-auto (Prive et al. 2020, 2023) estimates the point-normal prior's proportion of causal variants \(p\) and the heritability \(h^2\) within its Gibbs sampler, so nothing is tuned. One chain runs from each starting \(p\) in pInit. As in the LDpred2 vignette, a chain whose fitted correlations span less than 0.95 times the 95th percentile of the chains' spans is dropped, and the rest are averaged. A chain has diverged, and is dropped too, once a kept draw's squared effects sum past twice the squared marginal effects; bigsnpr's own sampler tests every iteration, burn-in included.

Usage

ldpred2AutoWeights(
  stat,
  ldReference,
  backend = NULL,
  h2Init = NULL,
  pInit = exp(seq(log(1e-04), log(0.2), length.out = 30)),
  nIter = 700,
  nBurnin = 500,
  shrinkCorr = 0.95,
  allowJumpSign = FALSE,
  pBounds = c(1e-05, 1),
  fitRetention = c("none", "slim", "full")
)

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 sampler reads z / sqrt(n) for each variant's own sample size.

ldReference

The LD reference the statistics' variants are matched to by stat$variantNames: a genotype panel (see readGenotypes), a GenotypeMixture or an LdMatrix.

backend

"bigsnpr" or "pecotmr". NULL (default) is pecotmr for a genotype reference and bigsnpr for an LdMatrix; pecotmr's sampler reads only a genotype reference.

h2Init

The starting heritability. NULL (default) estimates it for the window by LAVA, as estimateH2 does for QTL summary statistics, which reads the reference's sample size, so an LdMatrix must record nSamples. Floored at 0.001, which is also the start where LAVA gives no estimate.

pInit

The starting proportions of causal variants, one chain each.

nIter, nBurnin

Each chain's iterations, burn-in included, and its burn-in.

shrinkCorr

The LD shrinkage \(c\): the sampler reads \(cR + (1 - c)I\).

allowJumpSign

Keep a drawn effect whose sign differs from the current one; FALSE sets it to zero instead.

pBounds

The bounds on \(p\).

fitRetention

How much of the fit is kept as the "fit" attribute of the returned weights. "none" (default) attaches nothing. "slim" keeps chains, a table with one row per chain (its starting pInit, its estimates p and h2, whether it diverged and whether it was kept), and the kept chains' mean inclusion probabilities postp, p and h2. "full" also keeps beta, each chain's effects (variants x chains).

Value

A numeric vector of weights, the kept chains' mean effects. Unless fitRetention = "none", the fit is attached as attribute "fit".

Details

The defaults are the LDpred2 vignette's, without the MAF-dependent prior variance.

backend = "bigsnpr" runs bigsnpr's snp_ldpred2_auto() on the LD matrix, written to a temporary file for the call; for a genotype reference the matrix is formed from its genotypes. "pecotmr" runs pecotmr's port of the same sampler on a genotype reference's factor of the LD, never forming the matrix. bigsnpr's functions call foreach::registerDoSEQ(), which replaces any foreach parallel backend registered in the session. Both draw their seeds from R's random number generator, so set.seed() makes the weights reproducible.

Examples

if (requireNamespace("bigsnpr", quietly = TRUE)) {
  data(eqtlRegionExample)
  X <- eqtlRegionExample$X[, 1:30]
  y <- eqtlRegionExample$yRes
  stat <- list(
    cor = as.vector(cor(X, y)),
    n = rep(nrow(X), ncol(X)),
    variantNames = colnames(X)
  )
  ld <- LdMatrix(cor(X), nSamples = nrow(X))
  set.seed(1)
  w <- ldpred2AutoWeights(stat, ld, h2Init = 0.1, pInit = c(0.01, 0.1),
    fitRetention = "slim")
  attr(w, "fit")$chains
}
#> Loading required package: foreach
#> Loading required package: rngtools
#> # A tibble: 2 × 5
#>   pInit     p      h2 diverged kept 
#>   <dbl> <dbl>   <dbl> <lgl>    <lgl>
#> 1  0.01 0.342 0.00128 FALSE    TRUE 
#> 2  0.1  0.323 0.00153 FALSE    FALSE