Advertisement
Guest User

Untitled

a guest
Aug 24th, 2019
77
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.34 KB | None | 0 0
  1. library(pwr)
  2.  
  3. ratioSS <- function(n){
  4.  
  5. #calculate SESOI based on 33% power
  6. SESOI <- pwr.t.test(n = n, sig.level = 0.05,
  7. power = 0.33, type = "one.sample",
  8. alternative = "two.sided")$d
  9.  
  10. #calculate exact n needed to get 80% power on inferiority test
  11. n_needed <- pwr.t.test(d = SESOI, sig.level = 0.05,
  12. power = 0.8, type = "one.sample",
  13. alternative = "greater")$n
  14.  
  15. #use pnorm() to calculate power using 2.5x sample size heuristic
  16. nt <- 2.5 * n #use 2.5x heuristic
  17. #nt <- n_needed #check power when you used EXACTLY the n needed
  18. se <- 1/sqrt(nt)
  19. a <- qnorm(0.05, mean = SESOI, sd = se)
  20. trueMu <- 0
  21. tt <- pnorm(a, mean = 0, sd = se)
  22.  
  23. #return original experiment's sample size, exact n needed for 80% power, and the power using 2.5x heuristic
  24. return(list(oGss = n, n = n_needed/n, pow = tt))
  25.  
  26. }
  27.  
  28. #run for many samples
  29. ratios <- lapply(sort(rep(10:1000,100)), function(x) ratioSS(x))
  30.  
  31. #format output
  32. d <- do.call(rbind, ratios)
  33. d <- data.frame(d)
  34. names(d) <- c("OGss","n", "pow")
  35. d$OGss <- as.numeric(d$OGss)
  36. d$n <- as.numeric(d$n)
  37. d$pow <- as.numeric(d$pow)
  38.  
  39. #create Plots
  40. hist(d$pow)
  41. mean(d$pow)
  42.  
  43. plot(d$OGss, d$pow, type = "l", main = "Power using 2.5x vs sample size",
  44. xlab = "Original sample size", ylab = "power of infer. test using 2.5x")
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement