Statistical power calculation for one-sided exact binomial test

binom-test-power.R

##' Calculate the minimum number of samples required for a one-sided exact
##' binomial test to distinuish between two success probabilities
##' with specified alpha and power
##' 
##' @param alpha is desired false positive rate (probability of incorrectly rejecting the null)
##' @param required.power is the probability of correctly rejects the null when the alternate is true
##' @param p0 is the expected proportion of successes under the null
##' @param p1 is the proportion of successes under the alternate hypothesis
##' @return (Invisibly) A list containing the required sample size and the number of successful
##' trials required
##'
##' @references \itemize{
##' \item R mailing list 
##'     \url{http://r.789695.n4.nabble.com/Sample-size-calculations-for-one-sided-binomial-exact-test-td3964313.html}
##' \item Ahern (2001) \url{http://stat.ethz.ch/education/semesters/as2011/bio/ahernSampleSize.pdf}
##' \item Fleming (1982) \url{http://www.jstor.org/stable/2530297}
##' }
binom.test.power <- function(alpha=0.05, required.power=0.9, p0=0.9, p1=0.95) {
    ##' First calculate Fleming (1982) value (Biometrics)
    ##' This is an estimate of the required sample size
    fleming <- ((qnorm(required.power) * (p1 * (1 - p1)) ^ 1.5 +
                     qnorm(1 - alpha) * (p0 * (1 - p0)) ^ 0.5) / (p1 - p0)) ^ 2
    cat(sprintf("Approximate value by Fleming (1982) is %g\n", fleming))
    
    ##' As per Ahern (2001), search within the range 0.8 to 4 times
    ##' the Fleming value
    N.start <- floor(0.8 * fleming)
    N.stop <- ceiling(4 * fleming)
    cat(sprintf("Searching range %g to %g...\n", N.start, N.stop))
    N <- N.start:N.stop
    
    ## Required number of events at each sample size under null
    ## hypothesis (critical value)
    crit.val <- qbinom(p = 1 - alpha, size = N, prob = p0)
    
    ## Calculate beta (type II error) for each N under alternate
    ## hypothesis; 1 - beta is the power
    Power <- 1 - pbinom(crit.val, N, p1)
        
    ## Find the smallest sample size yielding at least the required power
    samp.size <- min(which(Power > required.power))
    
    ## Get and print the required number of events to reject the null
    ## given the sample size required
    cat(paste("Exact result is",crit.val[samp.size] + 1, "out of", N[samp.size], "\n"))
    invisible(list(successes=crit.val[samp.size] + 1, N=N[samp.size]))
}