NUll check for flashr2

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\)

########  spiked covariance model

set.seed(99)
T = 300
lambda_0 = rep(NA,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)
  X = X/sqrt(mean((X)^2))
  lambda_0[t] =  svd(X)$d[1]
  data =  flashr2::flash_set_data(X)
  g1= flashr2::flash_r1(data,verbose=F)
  if(sum(g1$EL) ==0 )  {flash2_rank_0[t] = 0}
}

set.seed(99)
lambda_1 = rep(NA,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, (2)/(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)
  X = X/sqrt(mean((X)^2))
  lambda_1[t] =  svd(X)$d[1]
  data = flashr2::flash_set_data(X)
  g1= flashr2::flash_r1(data,verbose=F)
  if(sum(g1$EL) ==0 )  {flash2_rank_1[t] = 0}
}



lambda_0_0 = lambda_0[which(flash2_rank_0 ==0)]
lambda_0_1 = lambda_0[which(flash2_rank_0 ==1)]


lambda_1_0 = lambda_1[which(flash2_rank_1 ==0)]
lambda_1_1 = lambda_1[which(flash2_rank_1 ==1)]

simple_roc <- function(labels, scores){
  labels <- labels[order(scores, decreasing=TRUE)]
  data.frame(TPR=cumsum(labels)/sum(labels), FPR=cumsum(!labels)/sum(!labels), labels)
}

lambda.roc = simple_roc(labels = c(rep(0,T),rep(1,T)),scores = c(lambda_0,lambda_1))
plot(lambda.roc$FPR,lambda.roc$TPR,type = "l",col = "black")
flash.roc = simple_roc(labels = c(rep(0,T),rep(1,T)),scores = c(flash2_rank_0,flash2_rank_1))
lines(flash.roc$FPR,flash.roc$TPR,type = "l",col = "red")

########  spiked covariance model

set.seed(99)
T = 300
lambda_0 = rep(NA,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)
  X = X/sqrt(mean((X)^2))
  lambda_0[t] =  svd(X)$d[1]
  data =  flashr2::flash_set_data(X)
  g1= flashr2::flash_r1(data,verbose=F)
  if(sum(g1$EL) ==0 )  {flash2_rank_0[t] = 0}
}

set.seed(99)
lambda_1 = rep(NA,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)
  X = X/sqrt(mean((X)^2))
  lambda_1[t] =  svd(X)$d[1]
  data = flashr2::flash_set_data(X)
  g1= flashr2::flash_r1(data,verbose=F)
  if(sum(g1$EL) ==0 )  {flash2_rank_1[t] = 0}
}



lambda_0_0 = lambda_0[which(flash2_rank_0 ==0)]
lambda_0_1 = lambda_0[which(flash2_rank_0 ==1)]


lambda_1_0 = lambda_1[which(flash2_rank_1 ==0)]
lambda_1_1 = lambda_1[which(flash2_rank_1 ==1)]

simple_roc <- function(labels, scores){
  labels <- labels[order(scores, decreasing=TRUE)]
  data.frame(TPR=cumsum(labels)/sum(labels), FPR=cumsum(!labels)/sum(!labels), labels)
}

lambda.roc = simple_roc(labels = c(rep(0,T),rep(1,T)),scores = c(lambda_0,lambda_1))
plot(lambda.roc$FPR,lambda.roc$TPR,type = "l",col = "black")
flash.roc = simple_roc(labels = c(rep(0,T),rep(1,T)),scores = c(flash2_rank_0,flash2_rank_1))
lines(flash.roc$FPR,flash.roc$TPR,type = "l",col = "red")

likelihood ratio

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] MASS_7.3-47     flashr2_0.2-3   workflowr_0.4.0 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     yaml_2.1.14       parallel_3.3.0   
[16] stringr_1.2.0     knitr_1.17        REBayes_0.85     
[19] rprojroot_1.2     grid_3.3.0        ggplot2_2.2.1    
[22] ashr_2.1-25       magrittr_1.5      backports_1.1.0  
[25] scales_0.4.1      codetools_0.2-15  htmltools_0.3.6  
[28] assertthat_0.2.0  softImpute_1.4    colorspace_1.3-2 
[31] labeling_0.3      stringi_1.1.5     Rmosek_7.1.2     
[34] lazyeval_0.2.0    munsell_0.4.3     doParallel_1.0.10
[37] pscl_1.4.9        truncnorm_1.0-7   SQUAREM_2016.8-2