Last updated: 2019-02-24

workflowr checks: (Click a bullet for more information)

  • ✔ R Markdown file: up-to-date

    Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

  • ✔ Environment: empty

    Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

  • ✔ Seed: set.seed(20190115)

    The command set.seed(20190115) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

  • ✔ Session information: recorded

    Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

  • ✔ Repository version: 2b87e47

    Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.

    Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated: 
    
Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    .sos/
    Ignored:    data/.DS_Store
    Ignored:    output/.DS_Store

Untracked files:
    Untracked:  docs/figure/test.Rmd/
    Untracked:  output/dscoutProblem475.rds
    Untracked:  output/dscoutProblem75.rds
    Untracked:  output/finemap_compare_random_data_null_dscout.rds
    Untracked:  output/finemap_compare_random_data_signal_dscout.rds
    Untracked:  output/finemap_compare_small_data_signal_dscout.rds
    Untracked:  output/finemap_compare_small_data_signal_dscout_RE8.rds
    Untracked:  output/random_data_100_sim_gaussian_null_1_get_sumstats_1_finemap_1.rds
    Untracked:  output/random_data_76.rds
    Untracked:  output/random_data_76_sim_gaussian_8.rds
    Untracked:  output/random_data_76_sim_gaussian_8_get_sumstats_1.rds
    Untracked:  output/small_data_42_sim_gaussian_36_get_sumstats_2_susie_z_2.rds
    Untracked:  output/small_data_92_sim_gaussian_30_get_sumstats_2_susie_z_2.rds

Unstaged changes:
    Modified:   analysis/SuSiEzRProblem42.Rmd
    Modified:   analysis/SuSiEzRProblem92.Rmd
    Modified:   analysis/SusieZPerformanceRE3.Rmd
    Modified:   output/dsc_susie_z_v_output.rds
    Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
Expand here to see past versions:
    File Version Author Date Message
    Rmd 2b87e47 zouyuxin 2019-02-24 wflow_publish(“analysis/SusieZPerformanceRE8.Rmd”)
    html ec8242a zouyuxin 2019-02-24 Build site.
    Rmd 588f363 zouyuxin 2019-02-24 wflow_publish(“analysis/SusieZPerformanceRE8.Rmd”)

DAP-z is not stable, so we exclude it from the simulation. It has two problems:

  1. The algorithm doesn’t stop.

  2. It uses more than 6GB memories.

There is no credible set for signals in FINEMAP. We generate the credible set of causal signals as the union of the variables included in the smallest set of causal configurations that already covered 95% of the total posterior probability.

We use datasets from dsc-finemap, which is from GTExV8 genotypes. The genotype matrix X is 574 by 1001.

Simulation with signals

We simulate a gaussian y under various number of causal variables, total percentage of variance explained (PVE) and whether the signals have equal effect. The reason I control the effect size is that if we random generate the effect size for the signals, some signals have large effect size by chance. Therefore these signals have larger PVE.

We fit SuSiE with L = 5, FINEMAP with max 5 causals.

library(dscrutils)
dscout = dscquery('output/finemap_compare_small_data_signal', target='method sim_gaussian.pve sim_gaussian.n_signal sim_gaussian.effect_weight score_susie.objective score_susie.converged score.total score.valid score.size score.signal_pip score_susie.purity score_susie.top score_susie.overlap ',group = c("score: score_susie score_finemap", "method: susie_z susie_z_init finemap"))

colnames(dscout) = c('DSC', 'method', 'output.file', 'pve', 'n_signal', 'effect_weight', 'score', 'objective', 'converged', 'total', 'valid', 'size', 'signal_pip', 'purity', 'top', 'overlap')
dscout$effect_weight[which(dscout$effect_weight == 'rep(1/n_signal, n_signal)')] = 'equal'
dscout$effect_weight[which(dscout$effect_weight != 'equal')] = 'notequal'
library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
library(knitr)
library(kableExtra)
library(susieR)
dscout = readRDS('output/finemap_compare_small_data_signal_dscout_RE8.rds')
dscout.susie = dscout[dscout$method == 'susie_z',]
dscout.susie.init = dscout[dscout$method == 'susie_z_init',]
dscout.finemap = dscout[dscout$method == 'finemap',]
  • Failed run:

There is one failed run for susie_z. The data has 5 signals and the total PVE is 0.8. Problem 42

  • Size of CS:
size.susie = aggregate(size~effect_weight+n_signal+pve, dscout.susie, FUN = function(x) round(mean(x), 2))
colnames(size.susie)[colnames(size.susie) == 'size'] <- 'size.susie'
size.susie.init = aggregate(size~effect_weight+n_signal+pve, dscout.susie.init, FUN = function(x) round(mean(x), 2))
colnames(size.susie.init)[colnames(size.susie.init) == 'size'] <- 'size.susie.init'
size.finemap = aggregate(size~effect_weight+n_signal+pve, dscout.finemap, FUN = function(x) round(mean(x), 2))
colnames(size.finemap)[colnames(size.finemap) == 'size'] <- 'size.finemap'

size = Reduce(function(...) merge(...),
       list(size.susie, size.susie.init, size.finemap))
size %>% kable() %>% kable_styling(bootstrap_options = c("striped", "condensed", "responsive"), full_width = F)
effect_weight n_signal pve size.susie size.susie.init size.finemap
equal 1 0.05 11.11 11.11 998.07
equal 1 0.10 4.22 4.23 992.40
equal 1 0.20 1.41 1.40 968.32
equal 1 0.60 0.27 0.27 942.25
equal 1 0.80 0.16 0.16 563.46
equal 10 0.05 11.02 10.75 996.23
equal 10 0.10 6.44 6.43 942.31
equal 10 0.20 2.32 2.34 681.95
equal 10 0.60 1.97 1.99 229.12
equal 10 0.80 1.67 1.71 147.27
equal 3 0.05 10.18 10.19 990.67
equal 3 0.10 4.72 4.79 905.46
equal 3 0.20 3.00 3.02 702.05
equal 3 0.60 0.90 0.91 427.31
equal 3 0.80 1.10 1.06 174.79
equal 5 0.05 9.08 9.32 992.99
equal 5 0.10 6.19 6.17 942.01
equal 5 0.20 2.64 2.64 580.85
equal 5 0.60 0.61 0.47 165.55
equal 5 0.80 0.88 0.82 99.75
notequal 1 0.05 11.11 11.11 998.08
notequal 1 0.10 4.22 4.23 992.31
notequal 1 0.20 1.41 1.40 967.49
notequal 1 0.60 0.27 0.27 941.80
notequal 1 0.80 0.16 0.16 565.49
notequal 10 0.05 8.38 8.29 992.35
notequal 10 0.10 4.16 4.15 967.02
notequal 10 0.20 2.51 2.52 902.51
notequal 10 0.60 0.45 0.45 746.36
notequal 10 0.80 0.28 0.27 664.63
notequal 3 0.05 8.46 8.47 993.11
notequal 3 0.10 6.02 5.93 972.49
notequal 3 0.20 5.41 5.41 885.79
notequal 3 0.60 0.49 0.49 712.27
notequal 3 0.80 0.16 0.16 548.23
notequal 5 0.05 5.96 5.96 994.36
notequal 5 0.10 2.80 2.80 979.18
notequal 5 0.20 2.53 2.53 927.49
notequal 5 0.60 0.73 0.73 855.17
notequal 5 0.80 0.48 0.48 644.88
  • Purity of CS:
purity.susie = aggregate(purity~effect_weight+n_signal+pve, dscout.susie, FUN = function(x) round(mean(x), 3))
colnames(purity.susie)[colnames(purity.susie) == 'purity'] <- 'purity.susie'
purity.susie.init = aggregate(purity~effect_weight+n_signal+pve, dscout.susie.init, FUN = function(x) round(mean(x), 3))
colnames(purity.susie.init)[colnames(purity.susie.init) == 'purity'] <- 'purity.susie.init'

purity = Reduce(function(...) merge(...),
       list(purity.susie, purity.susie.init))
purity %>% kable() %>% kable_styling(bootstrap_options = c("striped", "condensed", "responsive"), full_width = F)
effect_weight n_signal pve purity.susie purity.susie.init
equal 1 0.05 0.501 0.484
equal 1 0.10 0.422 0.421
equal 1 0.20 0.314 0.314
equal 1 0.60 0.138 0.138
equal 1 0.80 0.092 0.092
equal 10 0.05 0.231 0.215
equal 10 0.10 0.311 0.297
equal 10 0.20 0.252 0.257
equal 10 0.60 0.154 0.155
equal 10 0.80 0.133 0.133
equal 3 0.05 0.309 0.314
equal 3 0.10 0.338 0.337
equal 3 0.20 0.257 0.256
equal 3 0.60 0.146 0.146
equal 3 0.80 0.125 0.126
equal 5 0.05 0.245 0.244
equal 5 0.10 0.336 0.336
equal 5 0.20 0.260 0.260
equal 5 0.60 0.147 0.143
equal 5 0.80 0.132 0.132
notequal 1 0.05 0.501 0.484
notequal 1 0.10 0.422 0.421
notequal 1 0.20 0.314 0.314
notequal 1 0.60 0.138 0.138
notequal 1 0.80 0.092 0.092
notequal 10 0.05 0.460 0.459
notequal 10 0.10 0.427 0.428
notequal 10 0.20 0.307 0.306
notequal 10 0.60 0.178 0.178
notequal 10 0.80 0.132 0.133
notequal 3 0.05 0.411 0.408
notequal 3 0.10 0.393 0.394
notequal 3 0.20 0.308 0.308
notequal 3 0.60 0.111 0.111
notequal 3 0.80 0.084 0.084
notequal 5 0.05 0.512 0.512
notequal 5 0.10 0.482 0.482
notequal 5 0.20 0.376 0.376
notequal 5 0.60 0.188 0.188
notequal 5 0.80 0.124 0.124
  • Power:
valid = aggregate(valid ~ effect_weight + n_signal + pve, dscout.susie, sum)
total = aggregate(DSC~ effect_weight + n_signal + pve, dscout.susie, length)
total$total_true = total$DSC * total$n_signal
power.susie = merge(valid, total)
power.susie$power.susie = round(power.susie$valid/(power.susie$total_true), 3)
colnames(power.susie)[colnames(power.susie) == 'valid'] <- 'valid.susie'

valid = aggregate(valid ~ effect_weight + n_signal + pve, dscout.susie.init, sum)
total = aggregate(DSC~ effect_weight + n_signal + pve, dscout.susie.init, length)
total$total_true = total$DSC * total$n_signal
power.susie.init = merge(valid, total)
power.susie.init$power.susie.init = round(power.susie.init$valid/(power.susie.init$total_true), 3)
colnames(power.susie.init)[colnames(power.susie.init) == 'valid'] <- 'valid.susie.init'

valid = aggregate(valid ~ effect_weight + n_signal + pve, dscout.finemap, sum)
total = aggregate(DSC ~ effect_weight + n_signal + pve, dscout.finemap, length)
total$total_true = total$DSC * total$n_signal
power.finemap = merge(valid, total)
power.finemap$power.finemap = round(power.finemap$valid/(power.finemap$total_true),3)
colnames(power.finemap)[colnames(power.finemap) == 'valid'] <- 'valid.finemap'

power = Reduce(function(...) merge(...),
       list(power.susie, power.susie.init, power.finemap))
power = power[,-4]
power %>% kable() %>% kable_styling(bootstrap_options = c("striped", "condensed", "responsive","bordered"), full_width = F) %>% add_header_above(c(" ", " ", " "," ", "SuSiE z" = 2, "SuSiE z init" = 2,"FINEMAP" = 2)) %>% column_spec(c(6, 8, 10), bold = T)
SuSiE z
SuSiE z init
FINEMAP
effect_weight n_signal pve total_true valid.susie power.susie valid.susie.init power.susie.init valid.finemap power.finemap
equal 1 0.05 150 72 0.480 73 0.487 150 1.000
equal 1 0.10 150 65 0.433 65 0.433 150 1.000
equal 1 0.20 150 48 0.320 48 0.320 150 1.000
equal 1 0.60 150 21 0.140 21 0.140 139 0.927
equal 1 0.80 150 14 0.093 14 0.093 140 0.933
equal 10 0.05 1500 35 0.023 29 0.019 1494 0.996
equal 10 0.10 1500 66 0.044 64 0.043 1455 0.970
equal 10 0.20 1500 91 0.061 99 0.066 1320 0.880
equal 10 0.60 1500 69 0.046 131 0.087 1104 0.736
equal 10 0.80 1500 62 0.041 110 0.073 1004 0.669
equal 3 0.05 450 53 0.118 50 0.111 449 0.998
equal 3 0.10 450 98 0.218 95 0.211 446 0.991
equal 3 0.20 450 98 0.218 98 0.218 446 0.991
equal 3 0.60 450 61 0.136 64 0.142 438 0.973
equal 3 0.80 450 45 0.100 48 0.107 412 0.916
equal 5 0.05 750 49 0.065 50 0.067 748 0.997
equal 5 0.10 750 92 0.123 92 0.123 743 0.991
equal 5 0.20 750 123 0.164 123 0.164 715 0.953
equal 5 0.60 750 95 0.127 93 0.124 686 0.915
equal 5 0.80 750 84 0.112 86 0.115 653 0.871
notequal 1 0.05 150 72 0.480 73 0.487 150 1.000
notequal 1 0.10 150 65 0.433 65 0.433 150 1.000
notequal 1 0.20 150 48 0.320 48 0.320 150 1.000
notequal 1 0.60 150 21 0.140 21 0.140 138 0.920
notequal 1 0.80 150 14 0.093 14 0.093 142 0.947
notequal 10 0.05 1500 69 0.046 69 0.046 1492 0.995
notequal 10 0.10 1500 68 0.045 68 0.045 1471 0.981
notequal 10 0.20 1500 50 0.033 50 0.033 1414 0.943
notequal 10 0.60 1500 29 0.019 29 0.019 1323 0.882
notequal 10 0.80 1500 20 0.013 20 0.013 1246 0.831
notequal 3 0.05 450 63 0.140 62 0.138 449 0.998
notequal 3 0.10 450 62 0.138 61 0.136 449 0.998
notequal 3 0.20 450 57 0.127 57 0.127 443 0.984
notequal 3 0.60 450 17 0.038 17 0.038 428 0.951
notequal 3 0.80 450 13 0.029 13 0.029 400 0.889
notequal 5 0.05 750 81 0.108 81 0.108 748 0.997
notequal 5 0.10 750 78 0.104 78 0.104 743 0.991
notequal 5 0.20 750 66 0.088 66 0.088 732 0.976
notequal 5 0.60 750 32 0.043 32 0.043 716 0.955
notequal 5 0.80 750 20 0.027 20 0.027 667 0.889
  • FDR
valid = aggregate(valid ~ effect_weight + n_signal + pve, dscout.susie, sum)
total = aggregate(total~ effect_weight + n_signal + pve, dscout.susie, sum)
fdr.susie = merge(valid, total)
fdr.susie$fdr.susie = round((fdr.susie$total - fdr.susie$valid)/fdr.susie$total, 3)
colnames(fdr.susie)[colnames(fdr.susie) == 'valid'] <- 'valid.susie'
fdr.susie = fdr.susie[,-5]

valid = aggregate(valid ~ effect_weight + n_signal + pve, dscout.susie.init, sum)
total = aggregate(total~ effect_weight + n_signal + pve, dscout.susie.init, sum)
fdr.susie.init = merge(valid, total)
fdr.susie.init$fdr.susie.init = round((fdr.susie.init$total - fdr.susie.init$valid)/fdr.susie.init$total, 3)
colnames(fdr.susie.init)[colnames(fdr.susie.init) == 'valid'] <- 'valid.susie.init'
fdr.susie.init = fdr.susie.init[,-5]

valid = aggregate(valid ~ effect_weight + n_signal + pve, dscout.finemap, sum)
total = aggregate(size ~ effect_weight + n_signal + pve, dscout.finemap, sum)
fdr.finemap = merge(valid, total)
fdr.finemap$fdr.finemap = round((fdr.finemap$size - fdr.finemap$valid)/fdr.finemap$size, 3)
colnames(fdr.finemap)[colnames(fdr.finemap) == 'valid'] <- 'valid.finemap'
fdr.finemap = fdr.finemap[,-5]

fdr = Reduce(function(...) merge(...),
       list(fdr.susie, fdr.susie.init, fdr.finemap))

fdr %>% kable() %>% kable_styling(bootstrap_options = c("striped", "condensed", "responsive","bordered"), full_width = F) %>% add_header_above(c(" ", " ", " ", "SuSiE z" = 2, "SuSiE z init" = 2,"FINEMAP" = 2)) %>% column_spec(c(5, 7, 9), bold = T)
SuSiE z
SuSiE z init
FINEMAP
effect_weight n_signal pve valid.susie fdr.susie valid.susie.init fdr.susie.init valid.finemap fdr.finemap
equal 1 0.05 72 0.163 73 0.099 150 0.999
equal 1 0.10 65 0.097 65 0.044 150 0.999
equal 1 0.20 48 0.040 48 0.020 150 0.999
equal 1 0.60 21 0.000 21 0.000 139 0.999
equal 1 0.80 14 0.000 14 0.000 140 0.998
equal 10 0.05 35 0.327 29 0.341 1494 0.990
equal 10 0.10 66 0.241 64 0.256 1455 0.990
equal 10 0.20 91 0.278 99 0.277 1320 0.987
equal 10 0.60 69 0.179 131 0.115 1104 0.968
equal 10 0.80 62 0.151 110 0.091 1004 0.955
equal 3 0.05 53 0.361 50 0.375 449 0.997
equal 3 0.10 98 0.203 95 0.181 446 0.997
equal 3 0.20 98 0.084 98 0.084 446 0.996
equal 3 0.60 61 0.062 64 0.030 438 0.993
equal 3 0.80 45 0.062 48 0.040 412 0.984
equal 5 0.05 49 0.319 50 0.315 748 0.995
equal 5 0.10 92 0.264 92 0.264 743 0.995
equal 5 0.20 123 0.146 123 0.146 715 0.992
equal 5 0.60 95 0.021 93 0.021 686 0.972
equal 5 0.80 84 0.023 86 0.023 653 0.956
notequal 1 0.05 72 0.163 73 0.099 150 0.999
notequal 1 0.10 65 0.097 65 0.044 150 0.999
notequal 1 0.20 48 0.040 48 0.020 150 0.999
notequal 1 0.60 21 0.000 21 0.000 138 0.999
notequal 1 0.80 14 0.000 14 0.000 142 0.998
notequal 10 0.05 69 0.169 69 0.159 1492 0.990
notequal 10 0.10 68 0.081 68 0.093 1471 0.990
notequal 10 0.20 50 0.057 50 0.057 1414 0.990
notequal 10 0.60 29 0.000 29 0.000 1323 0.988
notequal 10 0.80 20 0.048 20 0.000 1246 0.988
notequal 3 0.05 63 0.203 62 0.215 449 0.997
notequal 3 0.10 62 0.127 61 0.129 449 0.997
notequal 3 0.20 57 0.066 57 0.066 443 0.997
notequal 3 0.60 17 0.056 17 0.056 428 0.996
notequal 3 0.80 13 0.188 13 0.188 400 0.995
notequal 5 0.05 81 0.147 81 0.147 748 0.995
notequal 5 0.10 78 0.082 78 0.082 743 0.995
notequal 5 0.20 66 0.043 66 0.043 732 0.995
notequal 5 0.60 32 0.030 32 0.030 716 0.994
notequal 5 0.80 20 0.130 20 0.130 667 0.993

Session information

sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14.3

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] susieR_0.6.4.0438 kableExtra_1.0.1  knitr_1.20        dplyr_0.7.8      

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.0        highr_0.7         compiler_3.5.1   
 [4] pillar_1.3.1      git2r_0.24.0      workflowr_1.1.1  
 [7] bindr_0.1.1       R.methodsS3_1.7.1 R.utils_2.7.0    
[10] tools_3.5.1       digest_0.6.18     lattice_0.20-38  
[13] evaluate_0.12     tibble_2.0.1      viridisLite_0.3.0
[16] pkgconfig_2.0.2   rlang_0.3.1       Matrix_1.2-15    
[19] rstudioapi_0.9.0  yaml_2.2.0        bindrcpp_0.2.2   
[22] stringr_1.3.1     httr_1.4.0        xml2_1.2.0       
[25] hms_0.4.2         grid_3.5.1        webshot_0.5.1    
[28] rprojroot_1.3-2   tidyselect_0.2.5  glue_1.3.0       
[31] R6_2.3.0          rmarkdown_1.11    purrr_0.2.5      
[34] readr_1.3.1       magrittr_1.5      whisker_0.3-2    
[37] backports_1.1.3   scales_1.0.0      htmltools_0.3.6  
[40] assertthat_0.2.0  rvest_0.3.2       colorspace_1.4-0 
[43] stringi_1.2.4     munsell_0.5.0     crayon_1.3.4     
[46] R.oo_1.22.0

This reproducible R Markdown analysis was created with workflowr 1.1.1