Multivariate Bayesian variable selection regression

Comparing MASH analysis with simple multivariate analysis

Previously we showed that even though univariate analysis with degenerated MASH model gives identical results to SuSiE as expected (with non-decreasing ELBO), in multivariate calculations the ELBO is not always non-decreasing. To investigate the issue we will 1) further simplify the problem and 2) isolate the problem to posterior calculation versus ELBO calculations and check which part is problematic. The best way to achieve both is to implement a simple Bayesian multivariate regression model with prior $b \sim MVN(0, U)$ where $U$ is known, instead of using MASH prior for $b$.

This feature is now implemented in the BayesianMultivariateRegression class with an interface added to the main function such that the code will be triggered when prior variance is a matrix.

With this feature (and with Yuxin's sharp eyes!!) we were able to identify an issue caused by inconsistent interface between mvsusieR::susie and mvsusieR::MashInitializer in handling residual variances. After patching the issue (interface fix still needs to be finalized) we are able to get consistent result between simple Bayesian multivariate regression and MASH; and MASH ELBO now increases.

Test the calculation agree with univariate code

In [1]:
library(mvsusieR)
set.seed(2)
L = 5

Loading required package: mashr
Loading required package: ashr
In [2]:
dat = mvsusie_sim1(r=1)

Run the simulated univariate data with SuSiE,

In [3]:
res = susieR::susie(dat$X,dat$y,L=L,scaled_prior_variance=0.2,estimate_residual_variance=F,estimate_prior_variance=F)
res$elbo
  1. -429.40495176513
  2. -409.549300389438
  3. -407.707739710942
  4. -407.7058409462
  5. -407.70583907902

Now run it with multivariate simple prior implementation in mvsusieR,

In [6]:
res = mvsusieR::mvsusie(dat$X,dat$y,L=L,prior_variance=0.2*var(dat$y),compute_objective=T,estimate_residual_variance=F,estimate_prior_variance=F)
dim(res$b2)

NULL

So it is confirmed that this implementation produces identical results as SuSiE runs.

Test multivariate analysis

In [35]:
set.seed(2)
dat = mvsusie_sim1(r=1)
In [42]:
devtools::load_all('~/GIT/software/mvsusieR')

Loading mvsusieR
In [43]:
res = mvsusie(dat$X,dat$y,L=L,prior_variance=dat$V,compute_objective=T,estimate_residual_variance=F,estimate_prior_variance=F)

[1]   1 500

Error in private$.pip[j] * matrix(private$.posterior_b2[j, ]) - tcrossprod(private$.pip[j] * : non-conformable arrays
Traceback:

1. mvsusie(dat$X, dat$y, L = L, prior_variance = dat$V, compute_objective = T, 
 .     estimate_residual_variance = F, estimate_prior_variance = F)
2. mvsusie_core(data, s_init, L, residual_variance, prior_variance, 
 .     prior_weights, estimate_residual_variance, estimate_prior_variance, 
 .     estimate_prior_method, precompute_covariances, compute_objective, 
 .     max_iter, tol, track_fit, verbose)   # at line 88-90 of file /home/gaow/GIT/software/mvsusieR/R/mvsusie.R
3. SuSiE_model$fit(data, prior_weights, estimate_prior_method, verbose)   # at line 244 of file /home/gaow/GIT/software/mvsusieR/R/mvsusie.R
4. private$SER[[l]]$compute_kl(d)   # at line 53 of file /home/gaow/GIT/software/mvsusieR/R/ibss_algorithm.R
5. private$compute_expected_loglik_partial(d)   # at line 25 of file /home/gaow/GIT/software/mvsusieR/R/single_effect_model.R
6. private$compute_expected_loglik_partial_multivariate(d)   # at line 38 of file /home/gaow/GIT/software/mvsusieR/R/single_effect_model.R
7. lapply(1:nrow(private$.posterior_b1), function(j) private$.pip[j] * 
 .     matrix(private$.posterior_b2[j, ]) - tcrossprod(private$.pip[j] * 
 .     private$.posterior_b1[j, ]))   # at line 54 of file /home/gaow/GIT/software/mvsusieR/R/single_effect_model.R
8. FUN(X[[i]], ...)

Here the ELBO is non-decreasing, as expected.

Compare with MASH based model

In [33]:
m_init = mvsusieR:::MashInitializer$new(list(dat$V), 1, prior_weight=1, null_weight=0)
In [34]:
res = mvsusieR::mvsusie(dat$X,dat$y,L=L,prior_variance=m_init,compute_objective=T,estimate_residual_variance=F,estimate_prior_variance=F)
dim(res$b2)

[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
[1] 500   2
[1]   2   2 500
  1. 5
  2. 500
  3. 2

The result agree with above. Now we use a different prior choice -- a diagonal prior covariance. We analyze with simple Bayesian multivariate regression,

In [9]:
res = mvsusieR::mvsusie(dat$X,dat$y,L=L,prior_variance=0.2*diag(ncol(dat$y)),compute_objective=T,estimate_residual_variance=F,estimate_prior_variance=F)
res$elbo
  1. -2174.05329261268
  2. -2166.18821093334
  3. -2161.71560740896
  4. -2158.70586771695
  5. -2156.96373258729
  6. -2156.8863212194
  7. -2156.88422194688
  8. -2156.8841094768

and with MASH based regression,

In [10]:
m_init = mvsusieR:::MashInitializer$new(list(0.2*diag(ncol(dat$y))), 1, prior_weight=1, null_weight=0,alpha=0)
In [11]:
res = mvsusieR::mvsusie(dat$X,dat$y,L=L,prior_variance=m_init,compute_objective=T,estimate_residual_variance=F,estimate_prior_variance=F)
res$elbo
  1. -2174.05329261268
  2. -2166.18821093334
  3. -2161.71560740896
  4. -2158.70586771696
  5. -2156.96373258729
  6. -2156.8863212194
  7. -2156.88422194688
  8. -2156.8841094768

So we are comfortable at this point that ELBO for multivariate analysis is done correctly.

Copyright © 2016-2020 Gao Wang et al at Stephens Lab, University of Chicago