Last updated: 2018-05-12

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(12345)

    The command set.seed(12345) 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: ddf9062

    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:    analysis/.DS_Store
        Ignored:    analysis/BH_robustness_cache/
        Ignored:    analysis/FDR_Null_cache/
        Ignored:    analysis/FDR_null_betahat_cache/
        Ignored:    analysis/Rmosek_cache/
        Ignored:    analysis/StepDown_cache/
        Ignored:    analysis/alternative2_cache/
        Ignored:    analysis/alternative_cache/
        Ignored:    analysis/ash_gd_cache/
        Ignored:    analysis/average_cor_gtex_2_cache/
        Ignored:    analysis/average_cor_gtex_cache/
        Ignored:    analysis/brca_cache/
        Ignored:    analysis/cash_deconv_cache/
        Ignored:    analysis/cash_fdr_1_cache/
        Ignored:    analysis/cash_fdr_2_cache/
        Ignored:    analysis/cash_fdr_3_cache/
        Ignored:    analysis/cash_fdr_4_cache/
        Ignored:    analysis/cash_fdr_5_cache/
        Ignored:    analysis/cash_fdr_6_cache/
        Ignored:    analysis/cash_plots_cache/
        Ignored:    analysis/cash_sim_1_cache/
        Ignored:    analysis/cash_sim_2_cache/
        Ignored:    analysis/cash_sim_3_cache/
        Ignored:    analysis/cash_sim_4_cache/
        Ignored:    analysis/cash_sim_5_cache/
        Ignored:    analysis/cash_sim_6_cache/
        Ignored:    analysis/cash_sim_7_cache/
        Ignored:    analysis/correlated_z_2_cache/
        Ignored:    analysis/correlated_z_3_cache/
        Ignored:    analysis/correlated_z_cache/
        Ignored:    analysis/create_null_cache/
        Ignored:    analysis/cutoff_null_cache/
        Ignored:    analysis/design_matrix_2_cache/
        Ignored:    analysis/design_matrix_cache/
        Ignored:    analysis/diagnostic_ash_cache/
        Ignored:    analysis/diagnostic_correlated_z_2_cache/
        Ignored:    analysis/diagnostic_correlated_z_3_cache/
        Ignored:    analysis/diagnostic_correlated_z_cache/
        Ignored:    analysis/diagnostic_plot_2_cache/
        Ignored:    analysis/diagnostic_plot_cache/
        Ignored:    analysis/efron_leukemia_cache/
        Ignored:    analysis/fitting_normal_cache/
        Ignored:    analysis/gaussian_derivatives_2_cache/
        Ignored:    analysis/gaussian_derivatives_3_cache/
        Ignored:    analysis/gaussian_derivatives_4_cache/
        Ignored:    analysis/gaussian_derivatives_5_cache/
        Ignored:    analysis/gaussian_derivatives_cache/
        Ignored:    analysis/gd-ash_cache/
        Ignored:    analysis/gd_delta_cache/
        Ignored:    analysis/gd_lik_2_cache/
        Ignored:    analysis/gd_lik_cache/
        Ignored:    analysis/gd_w_cache/
        Ignored:    analysis/knockoff_10_cache/
        Ignored:    analysis/knockoff_2_cache/
        Ignored:    analysis/knockoff_3_cache/
        Ignored:    analysis/knockoff_4_cache/
        Ignored:    analysis/knockoff_5_cache/
        Ignored:    analysis/knockoff_6_cache/
        Ignored:    analysis/knockoff_7_cache/
        Ignored:    analysis/knockoff_8_cache/
        Ignored:    analysis/knockoff_9_cache/
        Ignored:    analysis/knockoff_cache/
        Ignored:    analysis/knockoff_var_cache/
        Ignored:    analysis/marginal_z_alternative_cache/
        Ignored:    analysis/marginal_z_cache/
        Ignored:    analysis/mosek_reg_2_cache/
        Ignored:    analysis/mosek_reg_4_cache/
        Ignored:    analysis/mosek_reg_5_cache/
        Ignored:    analysis/mosek_reg_6_cache/
        Ignored:    analysis/mosek_reg_cache/
        Ignored:    analysis/pihat0_null_cache/
        Ignored:    analysis/plot_diagnostic_cache/
        Ignored:    analysis/poster_obayes17_cache/
        Ignored:    analysis/real_data_simulation_2_cache/
        Ignored:    analysis/real_data_simulation_3_cache/
        Ignored:    analysis/real_data_simulation_4_cache/
        Ignored:    analysis/real_data_simulation_5_cache/
        Ignored:    analysis/real_data_simulation_cache/
        Ignored:    analysis/rmosek_primal_dual_2_cache/
        Ignored:    analysis/rmosek_primal_dual_cache/
        Ignored:    analysis/seqgendiff_cache/
        Ignored:    analysis/simulated_correlated_null_2_cache/
        Ignored:    analysis/simulated_correlated_null_3_cache/
        Ignored:    analysis/simulated_correlated_null_cache/
        Ignored:    analysis/simulation_real_se_2_cache/
        Ignored:    analysis/simulation_real_se_cache/
        Ignored:    analysis/smemo_2_cache/
        Ignored:    data/LSI/
        Ignored:    docs/.DS_Store
        Ignored:    docs/figure/.DS_Store
        Ignored:    output/fig/
    
    Unstaged changes:
        Deleted:    analysis/cash_plots_fdp.Rmd
    
    
    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 cc0ab83 Lei Sun 2018-05-11 update
    html 0f36d99 LSun 2017-12-21 Build site.
    html 853a484 LSun 2017-11-07 Build site.
    html 1ea081a LSun 2017-07-03 sites
    html 899a5f8 LSun 2017-05-09 mosek improvements
    rmd 1df916a LSun 2017-05-09 mosek improvement
    html 27bcaf9 LSun 2017-05-08 writeups
    rmd 6bac69e LSun 2017-05-08 writeups

Introduction

We are experimenting different ways to improve the performance of Rmosek, using typical data sets of correlated \(z\) scores.

z <- read.table("../output/z_null_liver_777.txt")
sel <- c(32, 327, 355, 483, 778)
ord <- c(4, 9, 9, 4, 4)
source("../code/gdash.R")

Algorithm and variations

Recall that our biconvex optimization problem is as follows.

\[ \begin{array}{rl} \min\limits_{\pi,w} & -\sum\limits_{j = 1}^n\log\left(\sum\limits_{k = 1}^K\sum\limits_{l=1}^L\pi_k w_l f_{jkl} + \sum\limits_{k = 1}^K\pi_kf_{jk0}\right) - \sum\limits_{k = 1}^K\left(\lambda_k^\pi - 1\right)\log\pi_k + \sum\limits_{l = 1}^L\lambda_l^w\phi\left(\left|w_l\right|\right) \\ \text{subject to} & \sum\limits_{k = 1}^K\pi_k = 1\\ & \sum\limits_{l=1}^L w_l \varphi^{(l)}\left(z\right) + \varphi\left(z\right) \geq 0, \forall z\in \mathbb{R}\\ & w_l \text{ decay reasonably fast,} \end{array} \] where \(- \sum\limits_{k = 1}^K\left(\lambda_k^\pi - 1\right)\log\pi_k\) and \(+ \sum\limits_{l = 1}^L\lambda_l^w\phi\left(\left|w_l\right|\right)\) are to regularize \(\pi_k\) and \(w_l\), respectively.

This problem can be solved iteratively. Starting with an initial value, each step two convex optimization problems are solved.

With a given \(\hat w\), \(\hat\pi\) is solved by

\[ \begin{array}{rl} \min\limits_{\pi} & -\sum\limits_{j = 1}^n\log\left(\sum\limits_{k = 1}^K\pi_k \left(\sum\limits_{l=1}^L \hat w_l f_{jkl} + f_{jk0}\right)\right) - \sum\limits_{k = 1}^K\left(\lambda_k^\pi - 1\right)\log\pi_k\\ \text{subject to} & \sum\limits_{k = 1}^K\pi_k = 1 \ ,\\ \end{array} \]

which is readily available with functions in cvxr.

Meanwhile, with a given \(\hat\pi\), the optimization problem to obtain \(\hat w\) becomes

\[ \begin{array}{rl} \min\limits_{w} & -\sum\limits_{j = 1}^n\log\left(\sum\limits_{l=1}^Lw_l\left(\sum\limits_{k = 1}^K\hat\pi_k f_{jkl}\right) + \sum\limits_{k = 1}^K\hat\pi_kf_{jk0}\right) + \sum\limits_{l = 1}^L\lambda_l^w \phi \left(\left|w_l\right|\right) \\ \text{subject to} & \sum\limits_{l=1}^L w_l \varphi^{(l)}\left(z\right) + \varphi\left(z\right) \geq 0, \forall z\in \mathbb{R}\\ & w_l \text{ decay reasonably fast.} \end{array} \] The two constraints are important to prevent numerical instability. Yet they are not readily manifested as convex. Ideally, the regularization \(\sum\limits_{l = 1}^L\lambda_l^w\phi\left(\left|w_l\right|\right)\) would be able to capture the essence of these two constaints, and at the same time maintain the convexity. Different ideas will be explored.

In this part, we’ll mainly work on the “basic form” of the \(w\) optimization problem; that is, the optimization problem without regularization \(\sum\limits_{l = 1}^L\lambda_l^w\phi\left(\left|w_l\right|\right)\).

In the following, we’ll take a look at the basic form in its primal and dual formulations, optimized by w.mosek() and w.mosek.primal functions. First, they are applied to the correlated null, and their performance can be compared with the previous implementation by cvxr. Then, they are applied to the simulated non-null data sets, and the results are compared with those obtained by ASH, as well as the truth.

Basic form: primal

In its most basic form, the \(w\) estimation problem is as follows.

\[ \begin{array}{rl} \min\limits_{w} & -\sum\limits_{j = 1}^n \log\left(\sum\limits_{l=1}^Lw_l a_{jl} + a_{j0}\right) \\ \text{subject to} & \sum\limits_{l=1}^Lw_l a_{jl} + a_{j0} \geq 0, \forall j \ .\\ \end{array} \] Let the matrix \(A = \left[a_{jl}\right]\), the vector \(a = \left[a_{j0}\right]\), and we have its equivalent form,

\[ \begin{array}{rl} \min\limits_{w, g} & -\sum\limits_{j = 1}^n \log\left(g_j\right) \\ \text{subject to} & Aw + a = g \\ & g \geq 0\ . \end{array} \]

This problem can be coded by Rmosek as a “separable convex optimization” (SCOPT) problem.

Correlated null

Fitted w: -0.03653526 0.1999284 0.01047094 -0.02012681 
Time Cost in Seconds: 7.945 2.884 13.152 

Fitted w: 0.03394291 0.7345604 -0.1532514 0.1883912 -0.05504895 0.02297176 -0.006908406 0.001216017 -0.0003147927 
Time Cost in Seconds: 6.902 1.98 6.669 

Fitted w: 0.02262546 0.9213963 0.02180428 0.1760983 -0.01379674 0.003968132 -0.004904973 -0.0006085249 -0.0002990675 
Time Cost in Seconds: 5.375 1.932 6.116 

Fitted w: 0.04540484 -0.1275978 0.009150985 0.01004515 
Time Cost in Seconds: 19.294 2.066 14.959 

Fitted w: 0.005943196 0.3980709 -0.009222458 0.02630697 
Time Cost in Seconds: 6.636 1.954 6.815 

Signal \(+\) correlated error

Converged: FALSE 
Number of Iteration: 101 
Time Cost in Seconds: 2181.997 255.051 1644.499 

Converged: FALSE 
Number of Iteration: 101 
Time Cost in Seconds: 2634.559 296.365 1970.452 

Converged: TRUE 
Number of Iteration: 45 
Time Cost in Seconds: 553.188 127.588 540.733 

Converged: TRUE 
Number of Iteration: 49 
Time Cost in Seconds: 781.062 126.839 622.938 

Converged: FALSE 
Number of Iteration: 101 
Time Cost in Seconds: 2161.426 247.628 1565.216 

Basic form: dual

The primal \(w\) optimization problem has \(n + L\) variables and \(n\) constraints. When \(n\) is large, such as \(n = 10K\) in our simulations, the time cost is usually substantial. As the authors of REBayes pointed out, it’s better to work on the dual form, which is

\[ \begin{array}{rl} \min\limits_{v} & -\sum\limits_{j = 1}^n \log\left(v_j\right) + a^Tv \\ \text{subject to} & A^Tv = 0 \\ & v \geq 0\ . \end{array} \] The dual form has \(n\) variables and more importantly, only \(L\) constraints. As \(L \ll n\), the dual form is far more computationally efficient.

Rmosek provides optimal solutions to both primal and dual variables, so \(\hat w\) is readily available when working on \(v\) optimization. But we need to be careful on which dual variables to use, such as $sol$itr$suc, $sol$itr$slc, or others.

Correlated null

Fitted w: -0.03642797 0.1971637 0.01074601 -0.02181494 
Time Cost in Seconds: 0.485 0.049 0.388 

Fitted w: 0.03361434 0.7339189 -0.1544144 0.1875235 -0.0560049 0.02267003 -0.007151714 0.001187682 -0.0003321398 
Time Cost in Seconds: 0.524 0.009 0.372 

Fitted w: 0.02273596 0.9206369 0.02224234 0.1750273 -0.01339113 0.003607418 -0.004794795 -0.0006395078 -0.0002908422 
Time Cost in Seconds: 0.529 0.008 0.375 

Fitted w: 0.04544396 -0.1272823 0.008811108 0.009915311 
Time Cost in Seconds: 0.448 0.006 0.322 

Fitted w: 0.006084177 0.3976595 -0.009103231 0.02610567 
Time Cost in Seconds: 0.489 0.005 0.313 

Signal \(+\) correlated error

Converged: TRUE 
Number of Iteration: 21 
Time Cost in Seconds: 46.978 2.71 29.017 

Converged: TRUE 
Number of Iteration: 21 
Time Cost in Seconds: 45.406 5.343 52.01 

Converged: TRUE 
Number of Iteration: 16 
Time Cost in Seconds: 34.293 3.519 26.79 

Converged: TRUE 
Number of Iteration: 31 
Time Cost in Seconds: 64.198 3.902 39.897 

Converged: TRUE 
Number of Iteration: 12 
Time Cost in Seconds: 27.247 1.446 15.392 

Conclusion

Implementation by Rmosek can fully reproduce the results obtained by cvxr. Moreover, the dual form gives the same results as the primal dorm does, with only a fraction of time.

It seems hopeful that we’ll find a key to successfully tackle correlation in simultaneous inference, which has eluded Prof. Brad Efron for more than a decade.

Session information

sessionInfo()
R version 3.4.3 (2017-11-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.4

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/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] Rmosek_8.0.69     PolynomF_1.0-1    CVXR_0.95         REBayes_1.2      
[5] Matrix_1.2-12     SQUAREM_2017.10-1 EQL_1.0-0         ttutils_1.0-1    

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.16      knitr_1.20        whisker_0.3-2    
 [4] magrittr_1.5      workflowr_1.0.1   bit_1.1-12       
 [7] lattice_0.20-35   R6_2.2.2          stringr_1.3.0    
[10] tools_3.4.3       grid_3.4.3        R.oo_1.21.0      
[13] git2r_0.21.0      scs_1.1-1         htmltools_0.3.6  
[16] bit64_0.9-7       yaml_2.1.18       rprojroot_1.3-2  
[19] digest_0.6.15     gmp_0.5-13.1      ECOSolveR_0.4    
[22] R.utils_2.6.0     evaluate_0.10.1   rmarkdown_1.9    
[25] stringi_1.1.6     Rmpfr_0.6-1       compiler_3.4.3   
[28] backports_1.1.2   R.methodsS3_1.7.1



This reproducible R Markdown analysis was created with workflowr 1.0.1