Last updated: 2017-09-17

Code version: 78e3954

library(flashr2)
library(MASS)
########  spiked covariance model

set.seed(99)
T = 1000
lambda_0 = rep(NA,T)
flash_rank_0= rep(1,T)
flash2_rank_0= rep(1,T)
for(t in 1:T){
  P = 300
  N = 100
  Sigma_0 = diag(rep(1,P))
  X = mvrnorm(n = N,mu = rep(0,P), Sigma = Sigma_0)
  lambda_0[t] =  svd(X)$d[1]
  data = set_flash_data(X)
  g1= flash_r1(data,verbose=F)
  g2 = flashr::flash(X,objtype = "l",
                     ash_para_l = list(nullweight = 1 ),
                     ash_para_f = list(nullweight = 1 ))
  if(sum(g1$EL) ==0 )  {flash2_rank_0[t] = 0}
  if(sum(g2$l) ==0 )  {flash_rank_0[t] = 0}
}

set.seed(99)
T = 1000
lambda_1 = rep(NA,T)
flash_rank_1= rep(1,T)
flash2_rank_1= rep(1,T)
for(t in 1:T){
  P = 300
  N = 100
  Sigma_0 = diag(rep(1,P))
  nu = rnorm(P, 0, (6)/(sqrt(P)))
  index_nu = sample(seq(1:P),(P*0.95))
  nu[index_nu] = 0
  Sigma_1 = Sigma_0 + nu %*% t(nu)
  X = mvrnorm(n = N,mu = rep(0,P), Sigma = Sigma_1)
  lambda_1[t] =  svd(X)$d[1]
  data = set_flash_data(X)
  g1= flash_r1(data,verbose=F)
  g2 = flashr::flash(X,objtype = "l",
                     ash_para_l = list(nullweight = 1 ),
                     ash_para_f = list(nullweight = 1 ))
  if(sum(g1$EL) ==0 )  {flash2_rank_1[t] = 0}
  if(sum(g2$l) ==0 )  {flash_rank_1[t] = 0}
}


lambda_1_0_0 = lambda_0[which(flash_rank_0 ==0)]
lambda_1_0_1 = lambda_0[which(flash_rank_0 ==1)]
lambda_2_0_0 = lambda_0[which(flash2_rank_0 ==0)]
lambda_2_0_1 = lambda_0[which(flash2_rank_0 ==1)]

lambda_1_1_0 = lambda_1[which(flash_rank_1 ==0)]
lambda_1_1_1 = lambda_1[which(flash_rank_1 ==1)]
lambda_2_1_0 = lambda_1[which(flash2_rank_1 ==0)]
lambda_2_1_1 = lambda_1[which(flash2_rank_1 ==1)]

model

In this case the data sets are: \[H_0: X \sim N(0, \Sigma_0)\] \[H_1: X \sim N(0,\Sigma_1)\] where \(\Sigma_0 = I\) and \(\Sigma_1 = I + \nu \nu^T\)

library(ggplot2)
dat = data.frame(lambda = c(lambda_0,lambda_1), name = c(rep("null",length(lambda_0)), rep("alt", length(lambda_1)) ))
ggplot(dat,aes(x=lambda)) +
  geom_histogram(data=subset(dat,name == 'null'),fill = "red", alpha = 0.2) +
  geom_histogram(data=subset(dat,name == 'alt'),fill = "green", alpha = 0.2) 
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

  • red: null
  • green: alternative

flashr

library(ggplot2)
dat = data.frame(lambda = c(lambda_1_0_0,lambda_1_0_1,lambda_1_1_0,lambda_1_1_1), 
                 name = c(rep("null 0",length(lambda_1_0_0)),
                          rep("null 1",length(lambda_1_0_1)),
                          rep("alt 0",length(lambda_1_1_0)),
                          rep("alt 1",length(lambda_1_1_1))))
ggplot(dat,aes(x=lambda)) +
  geom_histogram(data=subset(dat,name == 'null 0'),fill = "red", alpha = 0.2) +
  geom_histogram(data=subset(dat,name == 'null 1'),fill = "black", alpha = 0.2) +
  geom_histogram(data=subset(dat,name == 'alt 0'),fill = "blue", alpha = 0.2) +
  geom_histogram(data=subset(dat,name == 'alt 1'),fill = "green", alpha = 0.2)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

  • red: null as rank 0
  • black: null as rank 1
  • blue: alternative as 0
  • green: alternative as 1

flashr2

library(ggplot2)
dat2 = data.frame(lambda = c(lambda_2_0_0,lambda_2_0_1,lambda_2_1_0,lambda_2_1_1), 
                 name = c(rep("null 0",length(lambda_2_0_0)),
                          rep("null 1",length(lambda_2_0_1)),
                          rep("alt 0",length(lambda_2_1_0)),
                          rep("alt 1",length(lambda_2_1_1))))
ggplot(dat,aes(x=lambda)) +
  geom_histogram(data=subset(dat2,name == 'null 0'),fill = "red", alpha = 0.2) +
  geom_histogram(data=subset(dat2,name == 'null 1'),fill = "black", alpha = 0.2) +
  geom_histogram(data=subset(dat2,name == 'alt 0'),fill = "blue", alpha = 0.2) +
  geom_histogram(data=subset(dat2,name == 'alt 1'),fill = "green", alpha = 0.2)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

  • red: null as rank 0
  • black: null as rank 1
  • blue: alternative as 0
  • green: alternative as 1

Session information

sessionInfo()
R version 3.3.0 (2016-05-03)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.12.6 (unknown)

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] ggplot2_2.2.1   MASS_7.3-47     flashr2_0.1-2   workflowr_0.4.0
[5] rmarkdown_1.6  

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.12      git2r_0.19.0      plyr_1.8.4       
 [4] iterators_1.0.8   tools_3.3.0       digest_0.6.12    
 [7] evaluate_0.10.1   tibble_1.3.3      gtable_0.2.0     
[10] lattice_0.20-35   rlang_0.1.2       Matrix_1.2-11    
[13] foreach_1.4.3     rstudioapi_0.6    yaml_2.1.14      
[16] parallel_3.3.0    stringr_1.2.0     knitr_1.17       
[19] REBayes_0.85      rprojroot_1.2     grid_3.3.0       
[22] irlba_2.2.1       flashr_0.1.1      ashr_2.1-25      
[25] magrittr_1.5      backports_1.1.0   scales_0.4.1     
[28] codetools_0.2-15  htmltools_0.3.6   assertthat_0.2.0 
[31] colorspace_1.3-2  labeling_0.3      stringi_1.1.5    
[34] Rmosek_7.1.2      lazyeval_0.2.0    munsell_0.4.3    
[37] doParallel_1.0.10 pscl_1.4.9        truncnorm_1.0-7  
[40] SQUAREM_2016.8-2

This R Markdown site was created with workflowr