smallTelescopesSESOI

library(pwr)

ratioSS <- function(n){

  #calculate SESOI based on 33% power
  SESOI <- pwr.t.test(n = n, sig.level = 0.05,
                      power = 0.33, type = "one.sample",
                      alternative = "two.sided")$d
  
  #calculate exact n needed to get 80% power on inferiority test
  n_needed <- pwr.t.test(d = SESOI, sig.level = 0.05,
                         power = 0.8, type = "one.sample",
                         alternative = "greater")$n
  
  #use pnorm() to calculate power using 2.5x sample size heuristic
  nt <- 2.5 * n #use 2.5x heuristic
  #nt <- n_needed #check power when you used EXACTLY the n needed
  se <- 1/sqrt(nt)
  a <- qnorm(0.05, mean = SESOI, sd = se)
  trueMu <- 0
  tt <- pnorm(a, mean = 0, sd = se)
  
  #return original experiment's sample size, exact n needed for 80% power, and the power using 2.5x heuristic
  return(list(oGss = n, n = n_needed/n, pow = tt))
  
}

#run for many samples
ratios <- lapply(sort(rep(10:1000,100)), function(x) ratioSS(x))

#format output
d <- do.call(rbind, ratios)
d <- data.frame(d)
names(d) <- c("OGss","n", "pow")
d$OGss <- as.numeric(d$OGss)
d$n <- as.numeric(d$n)
d$pow <- as.numeric(d$pow)

#create Plots
hist(d$pow)
mean(d$pow)

plot(d$OGss, d$pow, type = "l", main = "Power using 2.5x vs sample size",
     xlab = "Original sample size", ylab = "power of infer. test using 2.5x")