cashr comparison: ashr with uniform

Introduction

Investigation on some finer points:

1. Will the use of uniform instead of normal make the performance by ashr worse on FDP?

2. What’s happening when the median FDP given by cashr is better than other methods with \(\pi_0 = 0.99\).

ashr with uniform vs normal

source("../code/gdash_lik.R")
source("../code/count_to_summary.R")
library(ggplot2)

FDP <- function (FDR, qvalue, beta) {
  return(sum(qvalue <= FDR & beta == 0) / max(sum(qvalue <= FDR), 1))
}

TDP <- function (FDR, qvalue, beta) {
  return(sum(qvalue <= FDR & beta != 0) / sum(beta != 0))
}

boxplot.quantile <- function(x) {
  r <- quantile(x, probs = c(0.10, 0.25, 0.5, 0.75, 0.90))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  return(r)
}

boxplot.quantile.sq <- function (x) {
  r <- sqrt(quantile(x^2, probs = c(0.10, 0.25, 0.5, 0.75, 0.90)))
  names(r) <- c("ymin", "lower", "middle", "upper", "ymax")
  return(r)
}

mean.sq <- function (x) {
  r <- sqrt(mean(x^2))
  return(r)
}

mysqrt_trans <- function() {
  scales::trans_new("mysqrt", 
                    transform = base::sqrt,
                    inverse = function(x) ifelse(x<0, 0, x^2),
                    domain = c(0, Inf))
}

r <- readRDS("../data/liver.rds")

ngene <- 1e4
top_genes_index = function (g, X) {
  return(order(rowSums(X), decreasing = TRUE)[1 : g])
}
lcpm = function (r) {
  R = colSums(r)
  t(log2(((t(r) + 0.5) / (R + 1)) * 10^6))
}
Y = lcpm(r)
subset = top_genes_index(ngene, Y)
r = r[subset,]

nsamp <- 5
pi0 <- 0.9
q.vec <- seq(0.001, 0.20, by = 0.001)
q <- 0.1
method.name.FDR <- c("cashr", "BH", "qvalue", "ashr.n", "ashr.u", "locfdr")
method.col.FDR <- scales::hue_pal()(length(method.name.FDR))

sd.z <- sapply(z.list, sd)
Noise <- cut(sd.z, breaks = c(0, quantile(sd.z, probs = 1 : 2 / 3), Inf), labels = c("Pseudo Deflation", "In-between", "Pseudo Inflation"))
typical.noise <- order(sd.z)[floor(quantile(seq(sd.z), c(0.1, 0.5, 0.9)))]

##================================================================

FDP.list <- lapply(q.vec, function (q) {
  t(mapply(function(qvalue.mat, beta, q) {
    apply(qvalue.mat, 2, function (qvalue, q, beta) {
      FDP(q, qvalue, beta)
    }, q, beta)
  }, qvalue.list, beta.list, q))
})
names(FDP.list) <- q.vec

TDP.list <- lapply(q.vec, function(q) {
  t(mapply(function(qvalue.mat, beta, q) {
    apply(qvalue.mat, 2, function (qvalue, q, beta) {
      TDP(q, qvalue, beta)
    }, q, beta)
  }, qvalue.list, beta.list, q))
})
names(TDP.list) <- q.vec

z.list.sel <- z.list[typical.noise]
names(z.list.sel) <- c("Pseudo Deflation", "In-between", "Pseudo Inflation")
z.sep.ggdata <- reshape2::melt(z.list.sel, value.name = "z")
z.sep.ggdata$L1 <- factor(z.sep.ggdata$L1, levels = c("Pseudo Deflation", "In-between", "Pseudo Inflation"))
z.sep.plot <- ggplot(data = z.sep.ggdata, aes(x = z)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.2) +
  facet_wrap(~L1, nrow = 1) +
  stat_function(fun = dnorm, aes(color = "N(0,1)"), lwd = 1.5, show.legend = TRUE) +
  theme(axis.title.y = element_text(size = 15),
        axis.text.y = element_text(size = 10),
        strip.text = element_text(size = 15),
        legend.position = "left",
        legend.text = element_text(size = 12),
        legend.key = element_blank()
        )

FDP.q <- FDP.list[[which(q.vec == q)]]
FDP.q.noise.mat <- rbind.data.frame(
  cbind.data.frame(Noise = rep("All", length(Noise)),
                   FDP.q),
  cbind.data.frame(Noise,
                   FDP.q)
)
FDP.q.ggdata <- reshape2::melt(FDP.q.noise.mat, id.vars = c("Noise"), variable.name = "Method", value.name = "FDP")

################################################

TDP.q <- TDP.list[[which(q.vec == q)]]
TDP.q.noise.mat <- rbind.data.frame(
  cbind.data.frame(Noise = rep("All", length(Noise)),
                   TDP.q),
  cbind.data.frame(Noise,
                   TDP.q)
)
TDP.q.ggdata <- reshape2::melt(TDP.q.noise.mat, id.vars = c("Noise"), variable.name = "Method", value.name = "TDP")

FDP.sqrt.q.all.sep.plot <- ggplot(data = FDP.q.ggdata, aes(x = Method, y = FDP, fill = Method, color = Method)) +
  stat_summary(fun.data = boxplot.quantile.sq, geom = "boxplot", aes(width = 0.75), position = position_dodge(), show.legend = FALSE) +
  stat_summary(fun.y = mean.sq, geom = "point", position = position_dodge(width = 0.9), show.legend = FALSE, shape = 13, size = 3) +
  scale_x_discrete(limits = rev(levels(FDP.q.ggdata$Method))) +
  scale_y_continuous(trans = "mysqrt", breaks = c(0, 0.1, 0.2, 0.4, 0.6, 0.8)) +
  coord_flip() +
  scale_color_manual(labels = method.name.FDR, values = method.col.FDR) +
  scale_fill_manual(labels = method.name.FDR, values = alpha(method.col.FDR, 0.35)) +
  facet_wrap(~Noise, nrow = 1) +
  geom_hline(yintercept = q, col = "black", linetype = "dashed", size = 1) +
  labs(y = "FDP") +
  expand_limits(y = 0) +
  theme(plot.title = element_text(size = 12, hjust = 0),
        axis.title.y = element_blank(),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.text.x = element_text(size = 10),
        strip.text = element_text(size = 15),
        legend.position = "bottom",
        legend.background = element_rect(color = "grey"),
        legend.text = element_text(size = 12)
        )

Warning: Ignoring unknown aesthetics: width

TDP.q.all.sep.plot <- ggplot(data = TDP.q.ggdata, aes(x = Method, y = TDP, fill = Method, color = Method)) +
  stat_summary(fun.data = boxplot.quantile, geom = "boxplot", aes(width = 0.75), position = position_dodge(), show.legend = FALSE) +
  coord_flip() +
  stat_summary(fun.y = mean, geom = "point", position = position_dodge(width = 0.9), show.legend = FALSE, shape = 13, size = 3) +
  scale_x_discrete(limits = rev(levels(TDP.q.ggdata$Method))) +
  scale_color_manual(labels = method.name.FDR, values = method.col.FDR) +
  scale_fill_manual(labels = method.name.FDR, values = alpha(method.col.FDR, 0.35)) +
  facet_wrap(~Noise, nrow = 1) +
  labs(y = "TDP") +
  expand_limits(y = 0) +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_text(size = 15),
        axis.title.x = element_text(size = 15),
        axis.text.x = element_text(size = 10),
        strip.text = element_text(size = 15),
        legend.position = "bottom",
        legend.background = element_rect(color = "grey"),
        legend.text = element_text(size = 12)
        )

Warning: Ignoring unknown aesthetics: width

FDP.sqrt.TDP.q.sep.plot.save <- gridExtra::arrangeGrob(
  z.sep.plot +
    labs(title = "Typical Examples of Realized Correlated N(0,1) Noise") +
    scale_color_manual(values = "blue") +
    theme(legend.title = element_blank(),
         plot.margin = grid::unit(c(5.5, 5.5, 5.5, 67.5), "points"),
         plot.title = element_text(size = 12, hjust = 0.5),
         strip.text = element_text(size = 12),
         axis.title.x = element_blank(),
         axis.text.x = element_text(size = 10),
         axis.title.y = element_text(size = 12)
    ),
  gridExtra::arrangeGrob(
  FDP.sqrt.q.all.sep.plot +
    labs(title = "FDP") +
    theme(plot.title = element_text(size = 12, hjust = 0.5),
          axis.title.x = element_blank(),
          strip.text = element_text(size = 12),
          plot.margin = grid::unit(c(5.5, 5.5, 5.5, 5.5), "points"),
          axis.text.y = element_text(size = 12)
    ),
  TDP.q.all.sep.plot +
    labs(title = "TDP") +
    theme(plot.margin = grid::unit(c(5.5, 5.5, 5.5, 5.5), "points"),
          axis.title.x = element_blank(),
          strip.text = element_text(size = 12),
          plot.title = element_text(size = 12, hjust = 0.5),
          axis.text.y = element_text(size = 12)
    ),
  heights = c(1, 1),
  top = grid::textGrob(bquote(paste("Nominal FDR = ", .(q), " (", g[1], " is Gaussian; ", pi[0] == 0.9, ")")), gp = grid::gpar(fontsize = 12), hjust = 1.15)
  ),
  heights = c(1, 2)
)

ggsave("../output/fig/FDP.TDP.q.sep.unif.pdf", FDP.sqrt.TDP.q.sep.plot.save, height = 6, width = 8)

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

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     

other attached packages:
 [1] ggplot2_2.2.1     ashr_2.2-7        Rmosek_8.0.69    
 [4] PolynomF_1.0-2    CVXR_0.95         REBayes_1.3      
 [7] Matrix_1.2-14     SQUAREM_2017.10-1 EQL_1.0-0        
[10] ttutils_1.0-1    

loaded via a namespace (and not attached):
 [1] gmp_0.5-13.1      Rcpp_0.12.16      pillar_1.2.2     
 [4] plyr_1.8.4        compiler_3.4.3    git2r_0.21.0     
 [7] workflowr_1.1.1   R.methodsS3_1.7.1 R.utils_2.6.0    
[10] iterators_1.0.9   tools_3.4.3       digest_0.6.15    
[13] bit_1.1-13        tibble_1.4.2      gtable_0.2.0     
[16] evaluate_0.10.1   lattice_0.20-35   rlang_0.2.0      
[19] foreach_1.4.4     yaml_2.1.19       parallel_3.4.3   
[22] gridExtra_2.3     Rmpfr_0.7-0       ECOSolveR_0.4    
[25] stringr_1.3.1     knitr_1.20        rprojroot_1.3-2  
[28] bit64_0.9-7       grid_3.4.3        R6_2.2.2         
[31] rmarkdown_1.9     reshape2_1.4.3    magrittr_1.5     
[34] whisker_0.3-2     scales_0.5.0      MASS_7.3-50      
[37] backports_1.1.2   codetools_0.2-15  htmltools_0.3.6  
[40] scs_1.1-1         colorspace_1.3-2  labeling_0.3     
[43] stringi_1.2.2     lazyeval_0.2.1    munsell_0.4.3    
[46] pscl_1.5.2        doParallel_1.0.11 truncnorm_1.0-8  
[49] R.oo_1.22.0