Multivariate Bayesian variable selection regression

Single effect model sanity check

Check for agreement of different regression methods in VEM updates.

In [1]:
library(mvsusieR)
set.seed(2)
attach(mvsusie_sim1(r=1))

Loading required package: mashr
Loading required package: ashr
In [2]:
data = mvsusieR:::DenseData$new(X,y)
head(data$XtY)
A matrix: 6 × 1 of type dbl[,1]
-12.871105
45.080954
-32.900723
-41.838878
5.995901
27.611879
In [3]:
head(data$XtR)
A matrix: 6 × 1 of type dbl[,1]
-12.871105
45.080954
-32.900723
-41.838878
5.995901
27.611879

Univariate single effect regression

In [4]:
prior_var = 0.2 * as.numeric(var(y))
residual_var = as.numeric(var(y))
m1 = mvsusieR:::SingleEffectModel(mvsusieR:::BayesianSimpleRegression)$new(ncol(X), residual_var, prior_var)

Multivariate single effect regression

In [5]:
prior_covar = mvsusieR:::MashInitializer$new(list(0.2 * cov(y)), 1, alpha=0)
residual_covar = cov(y)
m2 = mvsusieR:::SingleEffectModel(mvsusieR:::MashRegression)$new(ncol(X), residual_covar, prior_covar)

Predictions after one fit

In [6]:
pred01 = m1$predict(data)
m1$fit(data)
pred1 = m1$predict(data)
In [7]:
pred02 = m2$predict(data)
m2$fit(data)
pred2 = m2$predict(data)
In [8]:
head(pred1)
A matrix: 6 × 1 of type dbl[,1]
-1.0398057
-1.6085843
-1.2046758
0.3079497
-2.3593036
-0.6760712
In [9]:
head(pred2)
A matrix: 6 × 1 of type dbl[,1]
-1.0398057
-1.6085843
-1.2046758
0.3079497
-2.3593036
-0.6760712
In [10]:
head(pred01)
A matrix: 6 × 1 of type dbl[,1]
0
0
0
0
0
0
In [11]:
head(pred02)
A matrix: 6 × 1 of type dbl[,1]
0
0
0
0
0
0
In [12]:
head(m1$pip)
  1. 4.21011575314195e-12
  2. 9.40264176321943e-12
  3. 6.24717097586617e-12
  4. 8.32853795168513e-12
  5. 3.98148063418841e-12
  6. 5.44318785610723e-12
In [13]:
head(m2$pip)
  1. 4.21011575314196e-12
  2. 9.40264176321945e-12
  3. 6.24717097586618e-12
  4. 8.32853795168514e-12
  5. 3.98148063418842e-12
  6. 5.44318785610724e-12
In [14]:
head(m1$posterior_b1 / m1$pip)
A matrix: 6 × 1 of type dbl[,1]
-0.06309365
0.22098507
-0.16127805
-0.20509254
0.02939167
0.13535235
In [15]:
head(m2$posterior_b1 / m2$pip)
A matrix: 6 × 1 of type dbl[,1]
-0.06309365
0.22098507
-0.16127805
-0.20509254
0.02939167
0.13535235
In [17]:
m1$lbf
18.1958934544426
In [18]:
m2$lbf
18.1958934544426
In [18]:
data$add_to_fitted(pred1)
data$compute_residual()
head(data$XtY)
-19.1055307
17.5706879
17.0295614
0.8703052
2.1676791
34.5852265
In [19]:
head(data$XtR)
0.7240101
8.8701269
9.5204835
-8.9370227
-1.3930653
28.8075105

Predictions after two fits

In [20]:
m1$fit(data)
pred1 = m1$predict(data)
In [21]:
m2$fit(data)
pred2 = m2$predict(data)
In [22]:
head(pred1)
0.008056959
-0.006026454
0.001472932
-0.001486319
-0.001913567
0.016121711
In [23]:
head(pred2)
0.008056959
-0.006026454
0.001472932
-0.001486319
-0.001913567
0.016121711
In [24]:
head(m1$pip)
  1. 0.00154586339736903
  2. 0.0016913068427532
  3. 0.00171474136974243
  4. 0.00169362641318872
  5. 0.0015483846584934
  6. 0.0040139083052572
In [25]:
head(m2$pip)
  1. 0.00154586339736903
  2. 0.0016913068427532
  3. 0.00171474136974243
  4. 0.00169362641318872
  5. 0.0015483846584934
  6. 0.0040139083052572
In [26]:
head(m1$posterior_b1)
5.486376e-06
7.353974e-05
8.002533e-05
-7.419597e-05
-1.057353e-05
5.668172e-04
In [27]:
head(m2$posterior_b1)
5.486376e-06
7.353974e-05
8.002533e-05
-7.419597e-05
-1.057353e-05
5.668172e-04
In [28]:
m1$lbf
-1.59617335538134
In [29]:
m2$lbf
-1.59617335538134
In [30]:
data$add_to_fitted(pred1)
data$compute_residual()
In [31]:
head(data$XtR)
0.6897485
8.8366233
9.3437291
-8.9747246
-1.4260285
28.4262408

Predictions after three fits

In [32]:
m1$fit(data)
pred1 = m1$predict(data)
In [33]:
m2$fit(data)
pred2 = m2$predict(data)
In [34]:
head(m1$posterior_b1)
5.268511e-06
7.380104e-05
7.886839e-05
-7.516685e-05
-1.091198e-05
5.498366e-04
In [35]:
head(m2$posterior_b1)
5.268511e-06
7.380104e-05
7.886839e-05
-7.516685e-05
-1.091198e-05
5.498366e-04
In [36]:
m1$lbf
-1.60418714571329
In [37]:
m2$lbf
-1.60418714571329

Check on main interface call

In [38]:
L = 10
In [39]:
residual_var = as.numeric(var(y))
scaled_prior_var = V[1,1] / residual_var
A = susie(X,y,L=L,prior_variance=scaled_prior_var,
          compute_objective=FALSE)
In [40]:
residual_cov = cov(y)
m_init = mvsusieR:::MashInitializer$new(list(V), 1, alpha = 0)
B = susie(X,y,L=L,prior_variance=m_init,compute_objective=FALSE)
In [44]:
max(abs(A$lbf - B$lbf))
1.4210854715202e-14
Copyright © 2016-2020 Gao Wang et al at Stephens Lab, University of Chicago