Multivariate Bayesian variable selection regression

Comparing multivariate MASH analysis using diagonal priors with univariate computations

This is to verify that the mvsusieR implementation is correct for the truthly multivariate computations. Previously I have only compared it for the degenerated case where the prior matrix is 1 by 1.

Simulate data

In [1]:
set.seed(1)
n = 1000
p = 1000
X = matrix(rnorm(n*p),nrow=n,ncol=p)
beta1 = beta2 = rep(0,p)
beta1[1:4] = runif(4)
beta2[1:4] = runif(4)
y1 = X %*% beta1 + rnorm(n)
y2 = X %*% beta2 + rnorm(n)
In [2]:
beta1[1:4]
  1. 0.614305987721309
  2. 0.320199134992436
  3. 0.205220706528053
  4. 0.831295673502609
In [3]:
beta2[1:4]
  1. 0.146828331751749
  2. 0.368926718831062
  3. 0.732725627720356
  4. 0.773408411303535

Run univariate computation

In [4]:
m1 = susieR::susie(X,y1,residual_variance = var(y1), estimate_residual_variance=F, scaled_prior_variance = 0.2, L=10)
bhat1 = coef(m1)
In [5]:
m2 = susieR::susie(X,y2,residual_variance = var(y2), estimate_residual_variance=F, scaled_prior_variance = 0.2, L=10)
bhat2 = coef(m2)
In [6]:
bhat1[2:5]
  1. 0.647417870851173
  2. 0.36269923511188
  3. 1.96971842692566e-11
  4. 0.81935068763822
In [7]:
bhat2[2:5]
  1. 3.72946361975251e-13
  2. 0.386342397093021
  3. 0.743487420194007
  4. 0.735604609996565

Run multivariate computation

In [8]:
prior_var = diag(0.2 * c(var(y1), var(y2)))
residual_var = diag(c(var(y1), var(y2)))
y = cbind(y1,y2)
data = mvsusieR:::DenseData$new(X,y)
In [9]:
prior_covar = mvsusieR:::MashInitializer$new(list(prior_var), 1, alpha=0)
In [10]:
m = mvsusieR::mvsusie(X, y, L=10, prior_variance=prior_covar)
In [11]:
dim(m$coef)
  1. 1001
  2. 2

Get lfsr

In [12]:
mvsusieR::mvsusie_get_cs_lfsr(m)
A matrix: 10 × 2 of type dbl[,2]
0.000000e+000.000000e+00
1.347289e-040.000000e+00
0.000000e+005.686698e-05
6.439294e-151.110223e-16
1.000000e+001.000000e+00
1.000000e+001.000000e+00
1.000000e+001.000000e+00
1.000000e+001.000000e+00
1.000000e+001.000000e+00
1.000000e+001.000000e+00

Get PIP per condition

In [13]:
pip_cond = mvsusieR:::mvsusie_get_pip_per_condition(m, prior_covar)
In [14]:
dim(pip_cond)
  1. 1000
  2. 2
In [15]:
head(pip_cond)
A matrix: 6 × 2 of type dbl[,2]
1.0000000001.000000000
1.0000000001.000000000
1.0000000001.000000000
1.0000000001.000000000
0.0035901180.003590118
0.0084418300.008441830
Copyright © 2016-2020 Gao Wang et al at Stephens Lab, University of Chicago