Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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")
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement