fdrate.R

#' False Discovery Rate (FDR)
#' see also http://www.sciencedirect.com/science/article/pii/S1053811901910377
#'
#'          Declared active  Declared inactive
#' Active       V_aa               V_ai       T_a
#' Inactive     V_ia               V_ii       T_i
#'              D_a                D_i         V
#'
#' instead of controlling the overall family wise error, we control
#' the proportion of false discoveries; that is, we want (V_aa) / (V_aa + V_ia)
#' to be smaller than q (which is most commonly picked as the alpha error 5%)
#'
#' @param: pvals: vector of p-values
#' @param: dist: string (default = 'dependent')
#' @param: q: numeric (maximum FDR a researcher is willing to risk)
#' @examples:
#' fdrte(c(0.02, 0.04, 0.03, 0.05))
fdrate <- function(pvals, dist = 'independent', q = 0.05) {
    V <- length(pvals)
    sorted <- sort(pvals)
    dist <- ifelse(dist == 'independent', 1, sum(1/1:length(pvals)))

    select <- sorted <= ((1:V/V) * q/dist)
    sorted[select]
}


#' Bonferroni correction
#'
#' @param: pvals: numeric vector
#' @param: alpha: numeric
#' bonferroni(c(.02, .04, .03, .05))
bonferroni <- function(pvals, alpha = .05) {
    select <- pvals <= (alpha / length(pvals))
    pvals[select]
}


#' FDR versus Bonferroni
#'
#' @param: d: numeric vector (Cohen's d)
#' @param: times: numeric (number of simulations)
#' @param: nr_tests: numeric (number of tests per simulation)
#' @param: n: numeric (sample size of each test)
#' @param: ...: additional arguments for fdrate
#' @examples:
#' sim(d = .5, times = 1000, nr_tests = 100, n = 50, q = .05)
sim <- function(d = 0.3, times = 100, nr_tests = 100, n = 50, ...) {
    res <- list('bonf' = 0, 'fdrate' = 0)

    for (i in 1:times) {
        pvals <- sapply(1:nr_tests, function(i) t.test(rnorm(n, d), rnorm(n))$p.value)
        res[['bonf']] <- res[['bonf']] + (length(bonferroni(pvals)) / nr_tests)
        res[['fdrate']] <- res[['fdrate']] + (length(fdrate(pvals, ...)) / nr_tests)
    }

    res
}

simres <- list('no_effect' = sim(d = 0),
               'small_effect' = sim(d = .1),
               'medium_effect' = sim(d = .3),
               'large_effect' = sim(d = .5))