Last updated: 2017-12-21

Code version: 6e42447

Introduction

Following previous simulation on randomly selected data sets, this time we first apply ash on \(z\) scores for every data set, and take a second look at those data sets on which ash makes most mistakes, or in other words, produces smallest \(\hat\pi_0\).

Simulation and ash

z = read.table("../output/z_null_liver_777.txt")
library(ashr)
n = dim(z)[1]
pihat0 = c()
for (i in 1:n) {
  ash.fit = ash(as.numeric(z[i, ]), 1, method = "fdr")
  pihat0[i] = get_pi0(ash.fit)
}

Histograms of \(z\) scores in 100 most “ash-hostile” data sets

I = order(pihat0)[1:100]
x = seq(-10, 10, 0.01)
y = dnorm(x)
k = 1
for (j in I) {
  cat("N0.", k, ": Data Set", j, "pihat0 =", pihat0[j])
  hist(as.numeric(z[j, ]), xlab = "z scores", freq = FALSE, ylim = c(0, 0.45), main = "10000 z scores, default")
  lines(x, y, col = "red")
  hist(as.numeric(z[j, ]), xlab = "z scores", freq = FALSE, ylim = c(0, 0.45), nclass = 100, main = "10000 z scores, 100 bins")
  lines(x, y, col = "red")
  k = k + 1
}
N0. 1 : Data Set 724 pihat0 = 0.01606004

N0. 2 : Data Set 33 pihat0 = 0.01680632

N0. 3 : Data Set 537 pihat0 = 0.01681555

N0. 4 : Data Set 522 pihat0 = 0.02083846

N0. 5 : Data Set 693 pihat0 = 0.02122591

N0. 6 : Data Set 840 pihat0 = 0.02316832

N0. 7 : Data Set 749 pihat0 = 0.02583276

N0. 8 : Data Set 269 pihat0 = 0.02787007

N0. 9 : Data Set 360 pihat0 = 0.0308234

N0. 10 : Data Set 23 pihat0 = 0.03139009

N0. 11 : Data Set 389 pihat0 = 0.03156173

N0. 12 : Data Set 923 pihat0 = 0.03239883

N0. 13 : Data Set 22 pihat0 = 0.03271534

N0. 14 : Data Set 885 pihat0 = 0.03413736

N0. 15 : Data Set 705 pihat0 = 0.03438913

N0. 16 : Data Set 40 pihat0 = 0.03668824

N0. 17 : Data Set 338 pihat0 = 0.0393062

N0. 18 : Data Set 379 pihat0 = 0.03958053

N0. 19 : Data Set 649 pihat0 = 0.04041784

N0. 20 : Data Set 627 pihat0 = 0.04041922

N0. 21 : Data Set 133 pihat0 = 0.04064386

N0. 22 : Data Set 915 pihat0 = 0.04083527

N0. 23 : Data Set 999 pihat0 = 0.04218153

N0. 24 : Data Set 182 pihat0 = 0.04223457

N0. 25 : Data Set 511 pihat0 = 0.04232901

N0. 26 : Data Set 122 pihat0 = 0.04301549

N0. 27 : Data Set 937 pihat0 = 0.04305909

N0. 28 : Data Set 771 pihat0 = 0.04316169

N0. 29 : Data Set 495 pihat0 = 0.04395344

N0. 30 : Data Set 41 pihat0 = 0.04491135

N0. 31 : Data Set 984 pihat0 = 0.04567195

N0. 32 : Data Set 3 pihat0 = 0.0465218

N0. 33 : Data Set 907 pihat0 = 0.04653847

N0. 34 : Data Set 355 pihat0 = 0.04750946

N0. 35 : Data Set 485 pihat0 = 0.04794909

N0. 36 : Data Set 56 pihat0 = 0.04829662

N0. 37 : Data Set 809 pihat0 = 0.04880257

N0. 38 : Data Set 404 pihat0 = 0.05059485

N0. 39 : Data Set 858 pihat0 = 0.05110069

N0. 40 : Data Set 77 pihat0 = 0.05151585

N0. 41 : Data Set 800 pihat0 = 0.05218407

N0. 42 : Data Set 562 pihat0 = 0.05247834

N0. 43 : Data Set 893 pihat0 = 0.05252202

N0. 44 : Data Set 403 pihat0 = 0.0531628

N0. 45 : Data Set 780 pihat0 = 0.05426134

N0. 46 : Data Set 538 pihat0 = 0.05493143

N0. 47 : Data Set 324 pihat0 = 0.05506971

N0. 48 : Data Set 176 pihat0 = 0.05513259

N0. 49 : Data Set 972 pihat0 = 0.05539611

N0. 50 : Data Set 644 pihat0 = 0.0557842

N0. 51 : Data Set 977 pihat0 = 0.0564768

N0. 52 : Data Set 467 pihat0 = 0.06112237

N0. 53 : Data Set 857 pihat0 = 0.06183555

N0. 54 : Data Set 783 pihat0 = 0.06208376

N0. 55 : Data Set 976 pihat0 = 0.06343862

N0. 56 : Data Set 968 pihat0 = 0.06385437

N0. 57 : Data Set 942 pihat0 = 0.06399678

N0. 58 : Data Set 402 pihat0 = 0.06468693

N0. 59 : Data Set 341 pihat0 = 0.06681535

N0. 60 : Data Set 411 pihat0 = 0.06812445

N0. 61 : Data Set 484 pihat0 = 0.06857855

N0. 62 : Data Set 476 pihat0 = 0.06944478

N0. 63 : Data Set 407 pihat0 = 0.07139176

N0. 64 : Data Set 871 pihat0 = 0.07236197

N0. 65 : Data Set 638 pihat0 = 0.07281058

N0. 66 : Data Set 232 pihat0 = 0.07317683

N0. 67 : Data Set 778 pihat0 = 0.07619716

N0. 68 : Data Set 349 pihat0 = 0.07627442

N0. 69 : Data Set 985 pihat0 = 0.07865805

N0. 70 : Data Set 259 pihat0 = 0.0794454

N0. 71 : Data Set 452 pihat0 = 0.079703

N0. 72 : Data Set 268 pihat0 = 0.08182634

N0. 73 : Data Set 526 pihat0 = 0.08340671

N0. 74 : Data Set 274 pihat0 = 0.08580661

N0. 75 : Data Set 685 pihat0 = 0.08592309

N0. 76 : Data Set 499 pihat0 = 0.08826948

N0. 77 : Data Set 982 pihat0 = 0.08854547

N0. 78 : Data Set 582 pihat0 = 0.08899861

N0. 79 : Data Set 953 pihat0 = 0.08900989

N0. 80 : Data Set 811 pihat0 = 0.08962902

N0. 81 : Data Set 575 pihat0 = 0.0946117

N0. 82 : Data Set 997 pihat0 = 0.09560788

N0. 83 : Data Set 231 pihat0 = 0.09686204

N0. 84 : Data Set 781 pihat0 = 0.09690435

N0. 85 : Data Set 853 pihat0 = 0.09889116

N0. 86 : Data Set 868 pihat0 = 0.09918104

N0. 87 : Data Set 502 pihat0 = 0.1001788

N0. 88 : Data Set 606 pihat0 = 0.1033893

N0. 89 : Data Set 912 pihat0 = 0.106048

N0. 90 : Data Set 659 pihat0 = 0.1065714

N0. 91 : Data Set 501 pihat0 = 0.1090034

N0. 92 : Data Set 296 pihat0 = 0.1115216

N0. 93 : Data Set 895 pihat0 = 0.1119629

N0. 94 : Data Set 153 pihat0 = 0.1148195

N0. 95 : Data Set 477 pihat0 = 0.1172838

N0. 96 : Data Set 769 pihat0 = 0.1187511

N0. 97 : Data Set 555 pihat0 = 0.1210902

N0. 98 : Data Set 475 pihat0 = 0.1273824

N0. 99 : Data Set 584 pihat0 = 0.1291777

N0. 100 : Data Set 949 pihat0 = 0.1308901

Take a closer look at the \(p\) values in 100 most “ash-hostile” data sets

We take a closer look at the behavior of \(p\) values in these most “ash-hostile” data sets, and also apply Benjamini-Hochberg’s FDR procedure on them, to see if BH is robust enough against the foes that make ash fail badly. False discoveries produced by BH are indicated as green points.

p = read.table("../output/p_null_liver_777.txt")
fd.bh = c()

for (i in 1 : n[1]) {
  p_BH = p.adjust(p[i, ], method = "BH")
  fd.bh[i] = sum(p_BH <= 0.05)
}
I = order(pihat0)[1:100]
x = seq(-10, 10, 0.01)
y = dnorm(x)
k = 1
m = dim(p)[2]
for (j in I) {
  cat("N0.", k, ": Data Set", j, "; pihat0 =", pihat0[j], "; Number of False Discoveries:", fd.bh[j])
  hist(as.numeric(z[j, ]), xlab = "z scores", freq = FALSE, ylim = c(0, 0.45), nclass = 100, main = "10000 z scores, 100 bins")
  lines(x, y, col = "red")
  qqnorm(as.numeric(z[j, ]), main = "Normal Q-Q plot for z scores")
  abline(0, 1)
  plot(sort(as.numeric(p[j, ])), xlab = "Order", ylab = "Ordered p value", main = "All p values")
  abline(0, 1 / m, col = "blue")
  abline(0, 0.05 / m, col = "red")
  legend("top", lty = 1, col = c("blue", "red"), c("Uniform", "BH, FDR = 0.05"))
  
  pj = sort(as.numeric(p[j, ]))

  plot(pj[1:max(100, fd.bh[j])], xlab = "Order", ylab = "Ordered p value", main = "Zoom-in to 100 smallest p values", ylim = c(0, pj[max(100, fd.bh[j])]))
  abline(0, 1 / m, col = "blue")
  points(pj[1:fd.bh[j]], col = "green", pch = 19)
  abline(0, 0.05 / m, col = "red")
  legend("top", lty = 1, col = c("blue", "red"), c("Uniform", "BH, FDR = 0.05"))
  
  plot(pj[pj <= 0.01], xlab = "Order", ylab = "Ordered p value", main = "Zoom-in to p values <= 0.01", ylim = c(0, 0.01))
  abline(0, 1 / m, col = "blue")
  points(pj[1:fd.bh[j]], col = "green", pch = 19)
  abline(0, 0.05 / m, col = "red")
  legend("top", lty = 1, col = c("blue", "red"), c("Uniform", "BH, FDR = 0.05"))
  
  plot(pj[pj <= 0.05], xlab = "Order", ylab = "Ordered p value", main = "Zoom-in to p values <= 0.05", ylim = c(0, 0.05))
  abline(0, 1 / m, col = "blue")
  points(pj[1:fd.bh[j]], col = "green", pch = 19)
  abline(0, 0.05 / m, col = "red")
  legend("top", lty = 1, col = c("blue", "red"), c("Uniform", "BH, FDR = 0.05"))
  
  k = k + 1
}
N0. 1 : Data Set 724 ; pihat0 = 0.01606004 ; Number of False Discoveries: 79

N0. 2 : Data Set 33 ; pihat0 = 0.01680632 ; Number of False Discoveries: 0

N0. 3 : Data Set 537 ; pihat0 = 0.01681555 ; Number of False Discoveries: 0

N0. 4 : Data Set 522 ; pihat0 = 0.02083846 ; Number of False Discoveries: 4

N0. 5 : Data Set 693 ; pihat0 = 0.02122591 ; Number of False Discoveries: 0

N0. 6 : Data Set 840 ; pihat0 = 0.02316832 ; Number of False Discoveries: 28

N0. 7 : Data Set 749 ; pihat0 = 0.02583276 ; Number of False Discoveries: 114

N0. 8 : Data Set 269 ; pihat0 = 0.02787007 ; Number of False Discoveries: 1

N0. 9 : Data Set 360 ; pihat0 = 0.0308234 ; Number of False Discoveries: 6

N0. 10 : Data Set 23 ; pihat0 = 0.03139009 ; Number of False Discoveries: 408

N0. 11 : Data Set 389 ; pihat0 = 0.03156173 ; Number of False Discoveries: 11

N0. 12 : Data Set 923 ; pihat0 = 0.03239883 ; Number of False Discoveries: 1

N0. 13 : Data Set 22 ; pihat0 = 0.03271534 ; Number of False Discoveries: 1

N0. 14 : Data Set 885 ; pihat0 = 0.03413736 ; Number of False Discoveries: 0

N0. 15 : Data Set 705 ; pihat0 = 0.03438913 ; Number of False Discoveries: 0

N0. 16 : Data Set 40 ; pihat0 = 0.03668824 ; Number of False Discoveries: 0

N0. 17 : Data Set 338 ; pihat0 = 0.0393062 ; Number of False Discoveries: 0

N0. 18 : Data Set 379 ; pihat0 = 0.03958053 ; Number of False Discoveries: 0

N0. 19 : Data Set 649 ; pihat0 = 0.04041784 ; Number of False Discoveries: 0

N0. 20 : Data Set 627 ; pihat0 = 0.04041922 ; Number of False Discoveries: 0

N0. 21 : Data Set 133 ; pihat0 = 0.04064386 ; Number of False Discoveries: 0

N0. 22 : Data Set 915 ; pihat0 = 0.04083527 ; Number of False Discoveries: 0

N0. 23 : Data Set 999 ; pihat0 = 0.04218153 ; Number of False Discoveries: 0

N0. 24 : Data Set 182 ; pihat0 = 0.04223457 ; Number of False Discoveries: 0

N0. 25 : Data Set 511 ; pihat0 = 0.04232901 ; Number of False Discoveries: 0

N0. 26 : Data Set 122 ; pihat0 = 0.04301549 ; Number of False Discoveries: 331

N0. 27 : Data Set 937 ; pihat0 = 0.04305909 ; Number of False Discoveries: 0

N0. 28 : Data Set 771 ; pihat0 = 0.04316169 ; Number of False Discoveries: 12

N0. 29 : Data Set 495 ; pihat0 = 0.04395344 ; Number of False Discoveries: 0

N0. 30 : Data Set 41 ; pihat0 = 0.04491135 ; Number of False Discoveries: 0

N0. 31 : Data Set 984 ; pihat0 = 0.04567195 ; Number of False Discoveries: 7

N0. 32 : Data Set 3 ; pihat0 = 0.0465218 ; Number of False Discoveries: 0

N0. 33 : Data Set 907 ; pihat0 = 0.04653847 ; Number of False Discoveries: 0

N0. 34 : Data Set 355 ; pihat0 = 0.04750946 ; Number of False Discoveries: 639

N0. 35 : Data Set 485 ; pihat0 = 0.04794909 ; Number of False Discoveries: 9

N0. 36 : Data Set 56 ; pihat0 = 0.04829662 ; Number of False Discoveries: 35

N0. 37 : Data Set 809 ; pihat0 = 0.04880257 ; Number of False Discoveries: 0

N0. 38 : Data Set 404 ; pihat0 = 0.05059485 ; Number of False Discoveries: 0

N0. 39 : Data Set 858 ; pihat0 = 0.05110069 ; Number of False Discoveries: 16

N0. 40 : Data Set 77 ; pihat0 = 0.05151585 ; Number of False Discoveries: 7

N0. 41 : Data Set 800 ; pihat0 = 0.05218407 ; Number of False Discoveries: 0

N0. 42 : Data Set 562 ; pihat0 = 0.05247834 ; Number of False Discoveries: 0

N0. 43 : Data Set 893 ; pihat0 = 0.05252202 ; Number of False Discoveries: 0

N0. 44 : Data Set 403 ; pihat0 = 0.0531628 ; Number of False Discoveries: 1

N0. 45 : Data Set 780 ; pihat0 = 0.05426134 ; Number of False Discoveries: 0

N0. 46 : Data Set 538 ; pihat0 = 0.05493143 ; Number of False Discoveries: 0

N0. 47 : Data Set 324 ; pihat0 = 0.05506971 ; Number of False Discoveries: 0

N0. 48 : Data Set 176 ; pihat0 = 0.05513259 ; Number of False Discoveries: 0

N0. 49 : Data Set 972 ; pihat0 = 0.05539611 ; Number of False Discoveries: 0

N0. 50 : Data Set 644 ; pihat0 = 0.0557842 ; Number of False Discoveries: 0

N0. 51 : Data Set 977 ; pihat0 = 0.0564768 ; Number of False Discoveries: 0

N0. 52 : Data Set 467 ; pihat0 = 0.06112237 ; Number of False Discoveries: 0

N0. 53 : Data Set 857 ; pihat0 = 0.06183555 ; Number of False Discoveries: 0

N0. 54 : Data Set 783 ; pihat0 = 0.06208376 ; Number of False Discoveries: 170

N0. 55 : Data Set 976 ; pihat0 = 0.06343862 ; Number of False Discoveries: 0

N0. 56 : Data Set 968 ; pihat0 = 0.06385437 ; Number of False Discoveries: 0

N0. 57 : Data Set 942 ; pihat0 = 0.06399678 ; Number of False Discoveries: 0

N0. 58 : Data Set 402 ; pihat0 = 0.06468693 ; Number of False Discoveries: 0

N0. 59 : Data Set 341 ; pihat0 = 0.06681535 ; Number of False Discoveries: 0

N0. 60 : Data Set 411 ; pihat0 = 0.06812445 ; Number of False Discoveries: 0

N0. 61 : Data Set 484 ; pihat0 = 0.06857855 ; Number of False Discoveries: 0

N0. 62 : Data Set 476 ; pihat0 = 0.06944478 ; Number of False Discoveries: 0

N0. 63 : Data Set 407 ; pihat0 = 0.07139176 ; Number of False Discoveries: 0

N0. 64 : Data Set 871 ; pihat0 = 0.07236197 ; Number of False Discoveries: 0

N0. 65 : Data Set 638 ; pihat0 = 0.07281058 ; Number of False Discoveries: 0

N0. 66 : Data Set 232 ; pihat0 = 0.07317683 ; Number of False Discoveries: 0

N0. 67 : Data Set 778 ; pihat0 = 0.07619716 ; Number of False Discoveries: 1

N0. 68 : Data Set 349 ; pihat0 = 0.07627442 ; Number of False Discoveries: 0

N0. 69 : Data Set 985 ; pihat0 = 0.07865805 ; Number of False Discoveries: 0

N0. 70 : Data Set 259 ; pihat0 = 0.0794454 ; Number of False Discoveries: 0

N0. 71 : Data Set 452 ; pihat0 = 0.079703 ; Number of False Discoveries: 0

N0. 72 : Data Set 268 ; pihat0 = 0.08182634 ; Number of False Discoveries: 0

N0. 73 : Data Set 526 ; pihat0 = 0.08340671 ; Number of False Discoveries: 0

N0. 74 : Data Set 274 ; pihat0 = 0.08580661 ; Number of False Discoveries: 2

N0. 75 : Data Set 685 ; pihat0 = 0.08592309 ; Number of False Discoveries: 0

N0. 76 : Data Set 499 ; pihat0 = 0.08826948 ; Number of False Discoveries: 0

N0. 77 : Data Set 982 ; pihat0 = 0.08854547 ; Number of False Discoveries: 0

N0. 78 : Data Set 582 ; pihat0 = 0.08899861 ; Number of False Discoveries: 0

N0. 79 : Data Set 953 ; pihat0 = 0.08900989 ; Number of False Discoveries: 0

N0. 80 : Data Set 811 ; pihat0 = 0.08962902 ; Number of False Discoveries: 0

N0. 81 : Data Set 575 ; pihat0 = 0.0946117 ; Number of False Discoveries: 1

N0. 82 : Data Set 997 ; pihat0 = 0.09560788 ; Number of False Discoveries: 1

N0. 83 : Data Set 231 ; pihat0 = 0.09686204 ; Number of False Discoveries: 0

N0. 84 : Data Set 781 ; pihat0 = 0.09690435 ; Number of False Discoveries: 0

N0. 85 : Data Set 853 ; pihat0 = 0.09889116 ; Number of False Discoveries: 0

N0. 86 : Data Set 868 ; pihat0 = 0.09918104 ; Number of False Discoveries: 0

N0. 87 : Data Set 502 ; pihat0 = 0.1001788 ; Number of False Discoveries: 0

N0. 88 : Data Set 606 ; pihat0 = 0.1033893 ; Number of False Discoveries: 0

N0. 89 : Data Set 912 ; pihat0 = 0.106048 ; Number of False Discoveries: 2

N0. 90 : Data Set 659 ; pihat0 = 0.1065714 ; Number of False Discoveries: 0

N0. 91 : Data Set 501 ; pihat0 = 0.1090034 ; Number of False Discoveries: 1

N0. 92 : Data Set 296 ; pihat0 = 0.1115216 ; Number of False Discoveries: 0

N0. 93 : Data Set 895 ; pihat0 = 0.1119629 ; Number of False Discoveries: 0

N0. 94 : Data Set 153 ; pihat0 = 0.1148195 ; Number of False Discoveries: 0

N0. 95 : Data Set 477 ; pihat0 = 0.1172838 ; Number of False Discoveries: 0

N0. 96 : Data Set 769 ; pihat0 = 0.1187511 ; Number of False Discoveries: 1

N0. 97 : Data Set 555 ; pihat0 = 0.1210902 ; Number of False Discoveries: 0

N0. 98 : Data Set 475 ; pihat0 = 0.1273824 ; Number of False Discoveries: 0

N0. 99 : Data Set 584 ; pihat0 = 0.1291777 ; Number of False Discoveries: 0

N0. 100 : Data Set 949 ; pihat0 = 0.1308901 ; Number of False Discoveries: 0

write.table(pihat0, file = "../output/pihat0_z_null_liver_777.txt", quote = FALSE, col.names = FALSE, row.names = FALSE)

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.2

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] compiler_3.4.3  backports_1.1.2 magrittr_1.5    rprojroot_1.3-1
 [5] tools_3.4.3     htmltools_0.3.6 yaml_2.1.16     Rcpp_0.12.14   
 [9] stringi_1.1.6   rmarkdown_1.8   knitr_1.17      git2r_0.20.0   
[13] stringr_1.2.0   digest_0.6.13   workflowr_0.8.0 evaluate_0.10.1

This R Markdown site was created with workflowr