simulate different probability distributions of p-values

pvalues-prob-sim.R

#----------------------------------------------------------------------------------------
# File: pvalues-prob-sim.R
# Author: Joseph Flanagan, reworking of script by Daniel Lakens
# email:   [email protected]
# Purpose: Function for demonstrating different probability distributions of p-values
#----------------------------------------------------------------------------------------

library(pwr)

sim_pvalues <- function(n, mean, sd, mu = 100, n_sims = 100000){
    p_values <- vapply(1:n_sims, FUN.VALUE = numeric(1), function(i){
        x <- rnorm(n = n, mean = mean, sd = sd)
        z <- t.test(x, mu = mu)
        p <- z$p.value
        return(p)
    })
    require(pwr)
    # determines M when power > 0. When power = 0, will set  M = 100.
    params <- pwr.t.test(d = (mean - 100)/sd, n = n, sig.level = 0.05, type = "one.sample", 
                        alternative = "two.sided")
    power <- params$power
    
    hist(p_values, breaks = 20, xlab = "p-values", ylab = "number of p-values\n", axes = FALSE,
        border = FALSE, 
        main = paste("p-value Distribution with", round(power * 100, digits = 1),"% Power"),
        col = "skyblue3",  xlim = c(0, 1), ylim = c(0, n_sims))
        axis(side = 1, at = seq(0, 1, 0.1), labels = seq(0, 1, 0.1))
        axis(side = 2, at = seq(0, n_sims, n_sims/4), labels = seq(0, n_sims, n_sims/4), las = 2)
        abline(h = n_sims/20, col = "red", lty = 3)
}

sim_pvalues(n = 26, mean = 106, sd = 15)