Calculate FDR using test statistics and null statistics, similar to `empPvals()` from qvalue package

getFDR.R

getFDR <- function(Stat, Stat0) {
  m <- length(Stat)
  n <- length(Stat0)
  v <- c(rep(T, m), rep(F, n))
  v <- v[order(c(Stat, Stat0), decreasing = T)]
  vSumT <- (m - cumsum(v == T)[v == T] + 1) / m
  vSumF <- (n - cumsum(v == F)[v == T]) / n
  fdr <- cummin(pmin(1, vSumF / vSumT))
  fdr <- fdr[rank(-Stat)]
  return(fdr)
}