bhack.R

library('BayesFactor')
# Bayesian version of http://rynesherman.com/phack.r
# see also http://osc.centerforopenscience.org/2014/07/02/phack/

bhack <- function(n = 30, hackrate = 10, m1 = 0, m2 = 0, sd1 = 1, sd2 = 1,
                  max_n = 100, bf_cutoff = 3, sims = 1000, graph = TRUE) {

    out <- matrix(NA, nrow = sims, ncol = 4)
    for (i in 1:sims) {
        hackcount <- 0
        group1 <- rnorm(n, m1, sd1)
        group2 <- rnorm(n, m2, sd2)
        initial_bf <- as.vector(ttestBF(group1, group2))[[1]]
        current_bf <- initial_bf

        while (current_bf < bf_cutoff & max_n > length(group1)) {
            hackcount <- hackcount + 1
            group1 <- c(group1, rnorm(hackrate, m1, sd1))
            group2 <- c(group2, rnorm(hackrate, m2, sd2))
            current_bf <- as.vector(ttestBF(group1, group2))[[1]]
        }

        n_added <- length(group1) - n
        out[i, ] <- c(initial_bf, hackcount, current_bf, n_added)
    }

    res <- data.frame(out)
    colnames(res) <- c("initial_bf", "hackcount", "final_bf", "n_added")
    initial_high_bf <- sum(res$initial_bf > bf_cutoff) / sims
    avg_hacks <- mean(res$hackcount)
    final_high_bf <- sum(res$final_bf > bf_cutoff) / sims
    avg_n_added <- mean(res$n_added)

    if(graph == TRUE) {
        op <- par(mfrow=c(2,1), las=1, font.main=1)
        init_sig <- res$initial_bf[res$initial_bf > bf_cutoff]
        final_sig <- res$final_bf[res$final_bf > bf_cutoff]
        #subtitle <- function(var) sprintf("Number of *significant* studies = %f (%f)", sims*var, var)

        hist(init_sig, xlab="BF-values", main="BF-curve for initial study",
             col="cyan", freq=FALSE); lines(density(init_sig))

        hist(final_sig, xlab="BF-values", main="BF-curve for hacked study",
             col="cyan", freq=FALSE); lines(density(final_sig))
    }

    cat("Proportion of Original Samples Statistically Significant =", initial_high_bf, "\n")
    cat("Proportion of Samples Statistically Significant After Hacking =", final_high_bf, "\n")
    cat("Average Number of Hacks Before Significant/Stopping =", avg_hacks, "\n")
    cat("Average N Added Before Significant/Stopping =", avg_n_added, "\n")
}