############################################################################################
 ##
## Stability of results from small RCT studies.
## Source: https://discourse.datamethods.org/t/stability-of-results-from-small-rcts/2287/11
## Original Author: R_cubed
## See also: ANDROMEDA-SHOCK trial discussion

quants <- c(0.00, 0.025, 0.05, 0.25, 0.50, 0.75, 0.95, 0.975, 1)

## 1. *p*-values under a true (point) null hypothesis.
##    In this scenario, *p*-values take on a uniform distribution. Calling this
##    function with any integer N is equivalent to simulating N studies with a
##    true difference of 0.
null_true_p <- function(n) {
  null   <- c(runif(n))
  result <- quantile(null, probs = quants)
  hist(null,
       freq = FALSE,
       main = expression(paste(H[0], " is True")),
       xlab = expression(paste(italic(p),"-values")))
  result
}


## 2. *p*-values under a true difference (standardized effect size of 0.53)
##    The median *p*-value varies based on how many studies are done, but is
##    frequently 0.75>p≤0.95 , with small numbers of studies (a common scenario
##    in meta-analysis). Note the skew the distribution. But skeptics will
##    almost always be able to find a large number of “contradictory” studies
##    going by “significance” alone.
null_false_p <- function(sample, sim_size)
{
  reject <- numeric(sim_size)
    for (i in 1:sim_size)
    {
      x <- rnorm(sample, 100, 15)
      y <- rnorm(sample, 108, 15)
      reject[i] <- t.test(x,y,alternative="greater")$p.value
      result    <- quantile(reject, probs = quants)
    }
    hist(reject,
         main  = expression(paste(H[alt], " is True")),
         freq  = FALSE,
         xlab  = expression(paste(italic(p),"-values")),
         sub   = "Effect Size=0.53, Sample Size=15")

    result
}

## 3. Meta-analysis under null
##    The following code simulates a meta-analysis based on *p*-value
##    combination using Stouffer’s method (ie. inverse normal transform).
##    95% of the (one-tailed) *t* statistics are within the range of
##    ±1.645; 95% are in the range of ±2.326, and 99.9% are in the range
##    of ±3.09 )
true_null_meta <- function(sample, sim_size)
{
  stouffer <- numeric(sim_size)
  for (i in 1:sim_size)
  {
    null        <- runif(sample)
    z_score     <- c(qnorm(null))
    stouffer[i] <- sum( z_score/sqrt(sample) )
    result      <- quantile(stouffer, probs = quants)
  }
  hist(stouffer,
       freq  = FALSE,
       main  = expression(paste(H[0], " is True")),
       xlab  = "Stouffer z-scores",
       sub   = paste("Number of Meta-Studies",sample))
  abline(v=qnorm(0.975),  col='red', lty=2)
  abline(v=-qnorm(0.975), col='red', lty=2)

  result
}


## 4. Meta-analysis with true effect
##    Finally, here is Stouffer’s method when the alternative is true, and
##    there is a difference. The magnitude of Stouffer’s statistic will depend
##    on the number of studies (as well as the magnitude of effect).
false_null_meta <- function(sample, sim_size) {
  z_scores <- numeric(sample)
  p_vals   <- numeric(sample)
  stouffer <- numeric(sim_size)
  for (i in 1:sim_size) {
    x <- rnorm(sample, 100, 15)
    y <- rnorm(sample, 108, 15)
    p_vals[i]   <- t.test(x,y, alternative = "greater")$p.value
    z_scores[i] <- c(qnorm(p_vals[i]))
    stouffer[i] <- sum( z_scores / sqrt(sample) )
    result_z    <- quantile(z_scores, probs = quants)
    result_s    <- quantile(stouffer, probs = quants)
  }
  hist(stouffer,
       freq  = FALSE,
       main  = expression(paste(H[alt], " is True")),
       xlab  = "Stouffer z-scores",
       sub   = paste("Number of Meta-Studies",sample))
  abline(v=qnorm(0.975),  col='red', lty=2)
  abline(v=-qnorm(0.975), col='red', lty=2)
  list(result_z=result_z, result_s=result_s)
}

N <- 1000
enter <- function() invisible(readline(prompt="Press [enter] to continue"))

cat("Distribution of *p*-values with effect size=0\n")
round(null_true_p(N), 3)
enter()

cat("Distribution of *p*-values with sample size 15 and effect size=0.53\n")
round(null_false_p(15, N), 3)
enter()

cat("Distribution of meta analysis z-scores with 15 studies\n")
round(true_null_meta(15, N), 3)
enter()

cat("Distribution of meta analysis z-scores with 15 studies and effect size=0.543\n")
false_null_meta(15, N)