Multivariate Bayesian variable selection regression

Visualization of M&M output

A prototype to plot various M&M output information.

Load a sample M&M output data-set

In [1]:
result_file = 'small_data_100_mixture01_1_oracle_generator_1_mnm_mixture01_1.rds'
meta_file = 'small_data_100_mixture01_1.pkl'
X_file = 'small_data_100.rds'
In [2]:
dat = readRDS(result_file)$result
names(dat)
  1. 'alpha'
  2. 'mu'
  3. 'mu2'
  4. 'KL'
  5. 'lbf'
  6. 'sigma2'
  7. 'V'
  8. 'elbo'
  9. 'niter'
  10. 'fitted'
  11. 'coef'
  12. 'null_index'
  13. 'mixture_weights'
  14. 'lfsr'
  15. 'intercept'
  16. 'sets'
  17. 'pip'
  18. 'm_init'
In [3]:
rownames(dat$coef) = NULL
In [4]:
meta = dscrutils::read_dsc(meta_file)
Y = meta$Y
X = readRDS(X_file)$X
In [5]:
univariate_res = lapply(1:ncol(Y), function(i) susieR:::univariate_regression(X,Y[,i]))
In [21]:
dat$bhat = do.call(cbind, lapply(1:ncol(Y), function(i) univariate_res[[i]]$betahat))
dat$shat = do.call(cbind, lapply(1:ncol(Y), function(i) univariate_res[[i]]$sebetahat))

SuSiE plot

We can directly call SuSiE plot function. Since it is simulated data we can add the true causal effects to it.

In [8]:
truth = meta$meta$true_coef
true_pos = as.integer(apply(truth, 1, sum) != 0)
true_idx = which(truth != 0, arr.ind = TRUE)
true_idx
rowcol
1491
3941
7861
1492
3942
7862
1493
3943
7863
1494
3944
7864
1495
3945
7865
In [9]:
truth[which(truth!=0,arr.ind=T)[,1],]
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
0.4118993 0.3934502 0.08806096 -0.072307930.05004165
In [10]:
dat$coef[-1,][true_idx]
  1. 5.54090525234181e-05
  2. 0.00897184235057198
  3. 0.0444660043879154
  4. 3.70363283541188e-05
  5. 0.00554368801446479
  6. 0.0307386813051094
  7. 1.26717037831625e-05
  8. 0.00200180257108648
  9. 0.00815300029442121
  10. -1.54978788368939e-05
  11. -0.00166610417021332
  12. -0.0074440800543281
  13. 7.12173907318122e-06
  14. 0.000838450772849928
  15. 0.00557786253834454
In [11]:
pdf('susie_plot_demo.pdf', width=10, height=5)
susieR::susie_plot(dat,y='PIP', main = 'Default SuSiE plot for cross-condition PIP', xlab = 'SNP positions', add_legend = T, b=true_pos)
dev.off()
png: 2
In [12]:
%preview susie_plot_demo.pdf -s png --dpi 100
> susie_plot_demo.pdf (11.0 KiB):

Bubble plot for estimated effect size

Effect size will be the color of the bubble:

In [13]:
p = mvsusieR::mvsusie_plot(dat)

Suggested PDF canvas width: 18 height: 4.5
In [14]:
pdf('bubble_demo.pdf', width = p$width, height = p$height)
print(p$plot)
dev.off()
png: 2
In [15]:
%preview bubble_demo.pdf -s png --dpi 100
> bubble_demo.pdf (6.4 KiB):
In [16]:
p = mvsusieR::mvsusie_plot(dat, cs_only=F)

Suggested PDF canvas width: 500.5 height: 4.5
In [17]:
pdf('bubble_demo_full.pdf', width = p$width, height = p$height)
print(p$plot)
dev.off()
png: 2
In [18]:
%preview bubble_demo_full.pdf -s png --dpi 100
> bubble_demo_full.pdf (31.8 KiB):

Bubble plot for original effect size

In [22]:
p = mvsusieR::mvsusie_plot(dat, original_sumstat = TRUE)

Suggested PDF canvas width: 18 height: 4.5
In [24]:
pdf('bubble_demo_z.pdf', width = p$width, height = p$height)
print(p$plot)
dev.off()
png: 2
In [25]:
%preview bubble_demo_z.pdf -s png --dpi 100
> bubble_demo_z.pdf (8.0 KiB):
In [28]:
p = mvsusieR::mvsusie_plot(dat, original_sumstat = TRUE, cs_only = FALSE)

Suggested PDF canvas width: 500.5 height: 4.5
In [30]:
pdf('bubble_demo_z_full.pdf', width = p$width, height = p$height)
print(p$plot)
dev.off()
png: 2
In [32]:
%preview bubble_demo_z_full.pdf -s png --dpi 100
> bubble_demo_z_full.pdf (87.4 KiB):
Copyright © 2016-2020 Gao Wang et al at Stephens Lab, University of Chicago