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 fe8f8f7 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 59fd661 LSun 2017-02-03 Build site.
    html 36c1e4c LSun 2017-02-03 Build site.
    Rmd d25a6e3 LSun 2017-02-03 step-down
    html d25a6e3 LSun 2017-02-03 step-down
    html c21d808 LSun 2017-02-02 Build site.
    Rmd 502f9e5 LSun 2017-02-02 Files commited by wflow_commit.
    html d575294 LSun 2017-01-31 Build site.
    Rmd 226434e LSun 2017-01-29 code
    html 226434e LSun 2017-01-29 code
    Rmd 6b6245e LSun 2017-01-29 url
    html 6b6245e LSun 2017-01-29 url
    Rmd ab2c030 LSun 2017-01-29 modify url
    html ab2c030 LSun 2017-01-29 modify url
    Rmd 1d187aa LSun 2017-01-29 First Commit
    html 1d187aa LSun 2017-01-29 First Commit

Last updated: 2018-05-15

Code version: 388e65e

Introduction

Suppose we have \(m + n\) observations of \((\hat\beta_j, \hat s_j)\) in two groups such that, with a pre-specified \(t\).

\[ \text{Group 1: }(\hat\beta_1, \hat s_1), \ldots, (\hat\beta_m, \hat s_m), \text{with } |\hat\beta_j/\hat s_j| \leq t \]

\[ \text{Group 2: }(\hat\beta_{m+1}, \hat s_{m+1}), \ldots, (\hat\beta_{m+n}, \hat s_{m+n}), \text{with } |\hat\beta_j/\hat s_j| > t \]

For Group 1, we’ll only use the information that for each one, \(|\hat\beta_j/\hat s_j| \leq t\); that is, they are moderate observations. For Group 2, we’ll use the full observation \((\hat\beta_j, \hat s_j)\).

Now we have usual ASH assumptions that for both groups,

\[ \begin{array}{c} \hat\beta_j | \hat s_j, \beta_j \sim N(\beta_j, \hat s_j^2)\\ \beta_j \sim \sum_k\pi_k N(0, \sigma_k^2) \end{array} \]

Now the estimate of \(\pi_k\) from the marginal probability of the data in both groups is

\[ \max_{\{\pi_k\}}\sum\limits_{j=1}^m\log(\sum_k\pi_k(2\Phi(\frac{t\hat s_j}{\sqrt{\sigma_k^2+\hat s_j^2}}) - 1)) + \sum\limits_{l=1}^{n}\log(\sum_k\pi_k N(\hat\beta_{m + l}; 0, \sigma_k^2 + \hat s_{m + l}^2)) \]

where \(\Phi(\cdot)\) is the cdf of \(N(0, 1)\), and \(N(\cdot; 0, \sigma_k^2 + \hat s_{m + l}^2)\) is the pdf of \(N(0, \sigma_k^2+\hat s_j^2)\).

Note that

  1. The second part of the objective function is the standard objective function for the regular ASH. The first part of the objective function is new, but it would be as if we are using a different likelihood, so it wouldn’t make any difference numerically.

  2. The first part of the objective,

\[ \max_{\{\pi_k\}}\sum\limits_{j=1}^m\log(\sum_k\pi_k(2\Phi(\frac{t\hat s_j}{\sqrt{\sigma_k^2+\hat s_j^2}}) - 1)) \]

and its special case, when we set \(\hat s_j \equiv 1\),

\[ m\log(\sum_k\pi_k(2\Phi(\frac{t}{\sqrt{\sigma_k^2+1}}) - 1)) \]

will move more weight to \(\hat\pi_1\). Customarily, we assume \(\sigma_1 \leq \cdots \leq \sigma_K\), so to maximize this first part alone will be essentially to maximize \(\hat\pi_1\). In other words, if we only have Group 1 and no “extreme” observations in Group 2, the estimated prior of \(\beta\) will be a point mass at \(0\). This is true no matter whether we are using \(z\)-scores (thus making \(\hat s_j \equiv 1\) effectively) or not, although using \(z\)-scores will make the case more obviously.

Implementation

Here is the main truncash function, utilizing autoselect.mixsd, mixIP, normalmix from ashr.

truncash

Simulation

The simplest case: \(g = 0.9\delta_0 + 0.1N(0, 1)\), independent observations

\(1000\) pairs of \((\hat\beta, \hat s \equiv 1)\) are simulated, so that \(900\) of them come from \(N(0, 1)\) and another 100 from \(N(\beta, 1)\) with \(\beta \sim N(0, 1)\). The \(1000\) pairs of observations are fed to both ash and truncash; \(\hat\pi_0\) and MSE are obtained. Results from 100 runs are plotted below.

source("../code/truncash.R")

set.seed(777)

sebetahat = rep(1, 1000)
t = qnorm(0.975)

pihat0_ash = pihat0_truncash = mse_ash = mse_truncash = mse_mle = lfsr0.05_ash = lfsr0.05_truncash = c()


for (i in 1:100) {
beta1 = rep(0, 900)
betahat1 = rnorm(length(beta1), beta1, 1)
beta2 = rnorm(100, 0, 1)
betahat2 = rnorm(length(beta2), beta2, 1)
beta = c(beta1, beta2)
betahat = c(betahat1, betahat2)

res.ash = ash(betahat, 1, mixcompdist = "normal", method = "fdr")
res.truncash = truncash(betahat, sebetahat, t)

pihat0_ash[i] = get_fitted_g(res.ash)$pi[1]
pihat0_truncash[i] = get_fitted_g(res.truncash)$pi[1]

mse_ash[i] = mean((get_pm(res.ash) - beta)^2)
mse_truncash[i] = mean((get_pm(res.truncash) - beta)^2)
mse_mle[i] = mean((betahat - beta)^2)

lfsr0.05_ash = mean(get_lfsr(res.ash) <= 0.05)
lfsr0.05_truncash = mean(get_lfsr(res.truncash) <= 0.05)
}

plot(pihat0_ash, pihat0_truncash, xlab = "ash", ylab = "truncash", main = expression(hat(pi)[0]),
     xlim = c(min(c(pihat0_ash, pihat0_truncash)), max(c(pihat0_ash, pihat0_truncash))),
     ylim = c(min(c(pihat0_ash, pihat0_truncash)), max(c(pihat0_ash, pihat0_truncash)))
     )
abline(0, 1, lty = 3, col = "blue")
abline(v = 0.9, lty = 3, col = "blue")
abline(h = 0.9, lty = 3, col = "blue")
points(0.9, 0.9, col = "red", cex = 1.5, pch = 19)

Expand here to see past versions of unnamed-chunk-1-1.png:
Version Author Date
0f36d99 LSun 2017-12-21
853a484 LSun 2017-11-07
d25a6e3 LSun 2017-02-03 
plot(abs(pihat0_ash - 0.9), abs(pihat0_truncash - 0.9), xlab = "ash", ylab = "truncash",
     main = expression(paste("|", hat(pi)[0] - 0.9, "|")), 
     xlim = c(0, max(c(abs(pihat0_ash - 0.9), abs(pihat0_truncash - 0.9)))),
     ylim = c(0, max(c(abs(pihat0_ash - 0.9), abs(pihat0_truncash - 0.9))))
     )
abline(0, 1, lty = 3, col = "blue")
points(0.9, 0.9, col = "red", cex = 1.5, pch = 19)

Expand here to see past versions of unnamed-chunk-1-2.png:
Version Author Date
0f36d99 LSun 2017-12-21
853a484 LSun 2017-11-07
d25a6e3 LSun 2017-02-03 
plot(mse_ash, mse_truncash, xlab = "ash", ylab = "truncash", main = "MSE",
     xlim = c(min(c(mse_ash, mse_truncash)), max(c(mse_ash, mse_truncash))),
     ylim = c(min(c(mse_ash, mse_truncash)), max(c(mse_ash, mse_truncash)))
     )
abline(0, 1, lty = 3, col = "blue")

Expand here to see past versions of unnamed-chunk-1-3.png:
Version Author Date
0f36d99 LSun 2017-12-21
853a484 LSun 2017-11-07
d25a6e3 LSun 2017-02-03

Session Information

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] SQUAREM_2017.10-1 ashr_2.2-2       

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   REBayes_1.2      
 [7] MASS_7.3-47       pscl_1.5.2        doParallel_1.0.11
[10] lattice_0.20-35   foreach_1.4.4     stringr_1.3.0    
[13] tools_3.4.3       parallel_3.4.3    grid_3.4.3       
[16] R.oo_1.21.0       git2r_0.21.0      htmltools_0.3.6  
[19] iterators_1.0.9   assertthat_0.2.0  yaml_2.1.18      
[22] rprojroot_1.3-2   digest_0.6.15     Matrix_1.2-12    
[25] codetools_0.2-15  R.utils_2.6.0     evaluate_0.10.1  
[28] rmarkdown_1.9     stringi_1.1.6     compiler_3.4.3   
[31] Rmosek_8.0.69     backports_1.1.2   R.methodsS3_1.7.1
[34] truncnorm_1.0-7

This reproducible R Markdown analysis was created with workflowr 1.0.1