Last updated: 2018-05-15

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: 388e65e

    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/
    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
    html e05bc83 LSun 2018-05-12 Update to 1.0
    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 b6bed27 LSun 2017-05-09 marginal null
    rmd 6264529 LSun 2017-05-09 marginal null
    html e5c73ae LSun 2017-05-05 marginal
    rmd 8206370 LSun 2017-05-05 marginal

This simulation can be seen as an enhanced version of a previous simulation.

Introduction

An assumption of using Gaussian derivatives to fit correlated null \(z\) scores is that each of these \(z\) scores should actually be null. That is, for \(n\) \(z\) scores \(z_1, \ldots, z_n\), although the correlation between \(z_i\) and \(z_j\) are not necessarily zero, the marginal distribution of \(z_i\), \(\forall i\), should be \(N\left(0, 1\right)\).

However, in practice, it’s not easy to check whether these correlated \(z\) scores are truly marginally \(N\left(0, 1\right)\). We’ve seen that their historgram could be far from normal. Further more, \(z\) scores in different data sets are distorted by different correlation structures. Therefore, we don’t have replicates here; that is, each data set is one single realization of a lot of random variables under correlation.

For our data sets in particular, let \(Z = \left[z_{ij}\right]_{m \times n}\) be the matrix of \(z\) scores. Each \(z_{ij}\) denotes the gene differential expression \(z\) score for gene \(j\) in the data set \(i\). Since all of these \(z\) scores are obtained from the same tissue, theoretically they should all be marginally \(N\left(0, 1\right)\).

Each row is a data set, consisting of \(10K\) realized \(z\) scores presumably marginally \(N\left(0, 1\right)\), whose empirical distribution distorted by correlation. If we plot the histogram of each row, it is grossly not \(N\left(0, 1\right)\) due to correlation. Therefore, it’s not easy to verify that they are truly marginally \(N\left(0, 1\right)\).

Here are two pieces of evidence that they are. Let’s take a look one by one, compared with the independent \(z\) scores case.

z.null <- read.table("../output/z_null_liver_777.txt")
n = ncol(z.null)
m = nrow(z.null)
set.seed(777)
z.sim = matrix(rnorm(m * n), nrow = m, ncol = n)

Row-wise: \(E\left[F_n\left(z\right)\right] = \Phi\left(z\right)\)

Let \(F_n^{R_i}\left(z\right)\) be the empirical CDF of \(p\) correlated \(z\) scores in row \(i\). For any \(i\), \(F_n^{R_i}\left(z\right)\) should be conspicuously different from \(\Phi\left(z\right)\), yet on average, \(E\left[F_n^{R_i}\left(z\right)\right]\) should be equal to \(\Phi\left(z\right)\), if all \(z\) scores are marginally \(N\left(0, 1\right)\).

In order to check that, we can borrow Prof. Michael Stein’s insight to look at the tail events, or empirical CDF.

For each row, let \(\alpha\) be a given probability level, \(z_\alpha = \Phi^{-1}\left(\alpha\right)\) be the associated quantile, and we record a number \(R_i^\alpha\) defined as follows.

If \(\alpha \leq 0.5\), \(R_\alpha^i\) is the number of \(z\) scores in row \(i\) that are smaller than \(z_\alpha\); otherwise, if \(\alpha > 0.5\), \(R_\alpha^i\) is the number of \(z\) scores in row \(i\) that are larger than \(z_\alpha\).

Defined this way, \(R_\alpha^i\) should be a sample from \(n \times F_n^{R_i}\left(z_\alpha\right)\) or \(n \times \left(1- F_n^{R_i}\left(z_\alpha\right)\right)\). We can check if \(E\left[F_n\left(z\right)\right] = \Phi\left(z\right)\) by looking at if the average \[ \bar R_\alpha \approx \begin{cases} n\Phi\left(z_\alpha\right) = n\alpha & \alpha \leq 0.5 \\ n\left(1-\Phi\left(z_\alpha\right)\right) = n\left(1 - \alpha\right) & \alpha > 0.5\end{cases} \ . \] We may also compare the frequencies of \(R_\alpha^i\) with their theoretical expected values \(m \times \text{Binomial}\left(n, \alpha\right)\) (in blue) assuming \(z_{ij}\) are independent.

Independent case: row-wise

Column-wise: closer to \(N\left(0, 1\right)\)

Each column of \(z\) should be seen as \(z\) scores of a non-differentially expressed gene in different data sets. Therefore, column-wise, the empirical distribution \(F_m^{C_j}\left(z\right)\) should be closer to \(\Phi\left(z\right)\) than \(F_m^{R_i}\left(z\right)\).

Similarly, we are plotting \(C_\alpha^i\), compared with their theoretical frequencies as follows.

Independent case: column wise

Conclusion

The empirical distribution and the indicated marginal distribution of the correlated null \(z\) scores are behaving not different from the expectation.

Row-wise, the number of tail observations averages to what would be expected from correlated marginally \(N\left(0, 1\right)\) random samples, validating Prof. Stein’s intuition.

Column-wise, the distribution of the number of tail observations is closer to normal, closer to what would be expected under corelated marginally \(N\left(0, 1\right)\). Moreover, the distribution seems unimodal, and peaked at \(m\alpha\) when \(\alpha \leq0.5\) or \(m\left(1-\alpha\right)\) when \(\alpha\geq0.5\). It suggests that the marginal distribution of the null \(z\) scores for a certain gene is usually centered at \(0\), and more often than not, close to \(N\left(0, 1\right)\).

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     

loaded via a namespace (and not attached):
 [1] workflowr_1.0.1   Rcpp_0.12.16      digest_0.6.15    
 [4] rprojroot_1.3-2   R.methodsS3_1.7.1 backports_1.1.2  
 [7] git2r_0.21.0      magrittr_1.5      evaluate_0.10.1  
[10] stringi_1.1.6     whisker_0.3-2     R.oo_1.21.0      
[13] R.utils_2.6.0     rmarkdown_1.9     tools_3.4.3      
[16] stringr_1.3.0     yaml_2.1.18       compiler_3.4.3   
[19] htmltools_0.3.6   knitr_1.20

This reproducible R Markdown analysis was created with workflowr 1.0.1