Last updated: 2017-12-21
Code version: 6e42447
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\).
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)
}
ash
-hostile” data setsI = 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
ash
-hostile” data setsWe 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)
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