Compute mvSuSiE-RSS TWAS weights from summary statistics

Multi-context summary-statistics analog of mvsusieWeights: extracts coefficients from an existing mvsusieR::mvsusie_rss fit, or fits one from stat$z (variants x conditions) and LD.

Usage

mvsusieRssWeights(
  stat,
  LD,
  mvsusieRssFit = NULL,
  priorVariance = NULL,
  residualVariance = NULL,
  L = 30,
  LGreedy = 5,
  retainFit = FALSE,
  ...
)

Arguments

stat

A list with z (matrix variants x conditions) and n (numeric vector or scalar).

LD

LD correlation matrix.

mvsusieRssFit

Optional pre-fitted mvsusieRss object.

priorVariance

Optional mvSuSiE prior variance specification. When NULL, mvsusieR::create_mixture_prior() is used with R = ncol(stat$z).

residualVariance

Optional residual covariance matrix.

L

Maximum number of single effects (default 30).

LGreedy

Initial greedy effect count (default 5).

retainFit

If TRUE, attaches the fitted object as an attribute.

...

Additional arguments forwarded to mvsusieR::mvsusie_rss.

Value

A numeric matrix of per-variant per-context weights (variants x conditions).

Details

Follows the *_rss_weights(stat, LD, ...) contract. Expects stat$z to be a numeric matrix (variants x conditions) and stat$n a per-context vector or scalar.

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(
  bhat = vapply(ss, `[`, numeric(1), 1L),
  shat = vapply(ss, `[`, numeric(1), 2L),
  z = vapply(ss, function(s) s[1] / s[2], numeric(1)),
  n = rep(nrow(X), ncol(X)))
LD <- cor(X)
fit <- fitMvsusieRss(Z = stat$z, R = LD, N = nrow(X),
  prior_variance = 1)
#> Eigendecomposition cache: K=1, common_cov=TRUE [mem: 0.90 GB]
#> Model initialized: J=30, R=1, L=10, K=1 [mem: 0.90 GB]
#> iter          ELBO       delta   sigma2      mem      V
#>    1     -588.3595           -   diag[1,1]   0.90 GB  [0 x 10]
#>    2     -588.3595    0.00e+00   diag[1,1]   0.90 GB  [0 x 10]  converged
mvsusieRssWeights(stat, LD, mvsusieRssFit = fit)
#>       [,1]
#>  [1,]    0
#>  [2,]    0
#>  [3,]    0
#>  [4,]    0
#>  [5,]    0
#>  [6,]    0
#>  [7,]    0
#>  [8,]    0
#>  [9,]    0
#> [10,]    0
#> [11,]    0
#> [12,]    0
#> [13,]    0
#> [14,]    0
#> [15,]    0
#> [16,]    0
#> [17,]    0
#> [18,]    0
#> [19,]    0
#> [20,]    0
#> [21,]    0
#> [22,]    0
#> [23,]    0
#> [24,]    0
#> [25,]    0
#> [26,]    0
#> [27,]    0
#> [28,]    0
#> [29,]    0
#> [30,]    0