plot_diagnostic in ashr

Introduction

This file documents some simulations done for a May 12, 2017, NHS presentation on the diagnostic plots for ASH implemented by ashr::plot_diagnostic().

library(ashr)
library(qvalue)

The concept of a Q-Q plot

set.seed(777)
x = rnorm(100)
qqnorm(x)
qqline(x)

y = runif(100)
qqnorm(y)
qqline(y)

Diagnostic plots for ASH applied to different data sets

beta = runif(10000)
betahat = beta + rnorm(10000)
fit = ash(betahat, 1, mixcompdist = "normal", method = "fdr")
plot_diagnostic(fit)

Press [enter] to see next plot

Press [enter] to see next plot

Press [enter] to see next plot

par(mfrow = c(2, 2))
fit = ash(rnorm(10000), 1, mixcompdist = "normal", method = "fdr")
plot_diagnostic(fit, breaks = 100)

Press [enter] to see next plot

Press [enter] to see next plot

Press [enter] to see next plot

beta = sample(c(rnorm(5000, 1.5, 1), rnorm(5000, -1.5, 1)))
betahat = beta + rnorm(10000)
fit = ash(betahat, 1, mixcompdist = "normal", method = "fdr")
plot_diagnostic(fit, breaks = 100)

Press [enter] to see next plot

Press [enter] to see next plot

Press [enter] to see next plot

The peculiarities of the correlated null

z = readRDS("../output/z_null_liver_777_select.RDS")
z = z$typical[[5]]

hist(z, breaks = 100, prob = TRUE, xlab = "Correlated N(0, 1) Z Scores", ylim = c(0, dnorm(0)), main = "Histogram of Correlated Z Scores")
lines(seq(-6, 6, 0.01), dnorm(seq(-6, 6, 0.01)), col = "red")

hist(z, breaks = 100, prob = TRUE, xlab = "Correlated N(0, 1) Z Scores", ylim = c(0, 0.1), main = "Histogram of Correlated Z Scores: Left Tail", xlim = c(-5, -2.5))
lines(seq(-6, 6, 0.01), dnorm(seq(-6, 6, 0.01)), col = "red")

pval = (1 - pnorm(abs(z))) * 2
q = qvalue(pval)
hist(pval, breaks = 100, main = "Histogram of Correlated Z Scores' p-values", xlab = "p-values", prob = TRUE)
abline(h = q$pi0, col = "red", lty = 2)
legend("topright", lty = 2, col = "red", expression(paste(hat(pi)[0], " by qvalue")))

Diagnostic plots for ASH applied to the correlated null.

library(ashr)
fit = ashr::ash(z, 1, mixcompdist = "normal", method = "fdr")
par(mfrow = c(2, 2))
plot_diagnostic(fit, breaks = 100)

Press [enter] to see next plot

Press [enter] to see next plot

Press [enter] to see next plot

Session information

sessionInfo()

R version 3.3.3 (2017-03-06)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS Sierra 10.12.4

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] qvalue_2.4.2 ashr_2.1-10 

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.10      knitr_1.15.1      magrittr_1.5     
 [4] REBayes_0.62      splines_3.3.3     MASS_7.3-45      
 [7] munsell_0.4.3     doParallel_1.0.10 pscl_1.4.9       
[10] colorspace_1.2-6  SQUAREM_2016.10-1 lattice_0.20-34  
[13] foreach_1.4.3     plyr_1.8.3        stringr_1.2.0    
[16] tools_3.3.3       parallel_3.3.3    grid_3.3.3       
[19] gtable_0.2.0      git2r_0.18.0      htmltools_0.3.5  
[22] iterators_1.0.8   assertthat_0.1    yaml_2.1.14      
[25] rprojroot_1.2     digest_0.6.11     Matrix_1.2-8     
[28] reshape2_1.4.1    ggplot2_2.1.0     codetools_0.2-15 
[31] evaluate_0.10     rmarkdown_1.3     stringi_1.1.2    
[34] Rmosek_7.1.3      scales_0.4.0      backports_1.0.5  
[37] truncnorm_1.0-7