seed_value <- 19710822 set.seed(seed_value) ## install.packages(c('pwr')) library(pwr) pwr.t.test(d = 45 / 30, power = 0.5, sig.level = 0.05, type = "two.sample", alternative = "two.sided") pwr.t.test(d = 45 / 30, power = 0.8) pwr.t.test(d = 45 / 30, power = 0.95) pwr.t.test(d = 13 / 30, power = 0.5) pwr.t.test(d = 13 / 30, power = 0.95) pwr.t.test(power = 0.8, n = 25) sample.sizes <- seq(5, 500, 5) effect.sizes <- numeric(length(sample.sizes)) for (i in 1:length(sample.sizes)) { effect.sizes[i] = pwr.t.test(n = sample.sizes[i], power = 0.5)$d } plot(sample.sizes, effect.sizes, type = "l") sample.sizes <- seq(5, 500, 5) effect.sizes <- numeric(length(sample.sizes)) for (i in 1:length(sample.sizes)) { effect.sizes[i] = pwr.t.test(n = sample.sizes[i], power = 0.5)$d } plot(sample.sizes, effect.sizes, type = "l") for (i in 1:length(sample.sizes)) { effect.sizes[i] = pwr.t.test(n = sample.sizes[i], power = 0.6)$d } lines(sample.sizes, effect.sizes, col='red') power.plot.1 <- function() { sample.sizes <- seq(5, 500, 10) powers = c(0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99) effect.sizes <- array(numeric(), c(length(sample.sizes), length(powers))) for (i in 1:length(sample.sizes)) { for (j in 1:length(powers)) { effect.sizes[i, j] = pwr.t.test(n = sample.sizes[i], power = powers[j])$d } } xrange <- c(floor(min(sample.sizes)), ceiling(max(sample.sizes))) yrange <- c(floor(min(effect.sizes)), ceiling(max(effect.sizes))) plot(xrange, yrange, xlab = 'Sample Size', ylab = 'Effect Size', type = 'n') for (j in 1:length(powers)) { lines(sample.sizes, effect.sizes[ ,j], col = j, lwd = 2) } legend('topright', title = 'Power', as.character(powers), lwd = 3, col = 1:j) } power.plot.1() power.plot.2 <- function() { sample.sizes <- c(3:15, seq(20, 50, 5)) effect.sizes <- seq(0.6, 1.6, 0.2) precisions <- effect.sizes * 35 powers <- array(numeric(), c(length(sample.sizes), length(effect.sizes))) for (i in 1:length(sample.sizes)) { for (j in 1:length(effect.sizes)) { powers[i, j] = pwr.t.test(n = sample.sizes[i], d = effect.sizes[j])$power } } xrange <- c(floor(min(sample.sizes)), ceiling(max(sample.sizes))) yrange <- c(0, 1) plot(xrange, yrange, xlab = 'Sample Size (each group)', ylab = 'Power', type = 'n') for (j in 1:length(effect.sizes)) { lines(sample.sizes, powers[ ,j], col = j, lwd = 3) } legend('bottomright', title = 'Precision', paste(as.character(precisions), "ohm*cm^2"), lwd = 3, col = 1:j) } power.plot.2() pwr.t.test(power = 0.80, d = 0.5)$n pwr.t2n.test(n1 = 48, power = 0.80, d = 0.5)$n2 pwr.t2n.test(n1 = 30, power = 0.80, d = 0.5)$n2 pwr.t.test.sim <- function(d = 1, n = 100, sig.level = 0.05, trial.count = 10000) { trials <- 1:trial.count p.values <- numeric(trial.count) for (i in trials) { x <- rnorm(n) y <- rnorm(n, mean = d) p.values[i] <- t.test(x, y)$p.value } hist(p.values, breaks = seq(0, 1, sig.level)) power <- sum(p.values <= sig.level) / trial.count return(power) } pwr.t.test(n = 50, d = 0.5) pwr.t.test.sim(n = 50, d = 0.5) pwr.t2sd.test.sim <- function(n = 100, mean1 = 0, sd1 = 1.0, mean2 = 1, sd2 = 1.0, sig.level = 0.05, trial.count = 10000) { trials <- 1:trial.count p.values <- numeric(trial.count) for (i in trials) { x <- rnorm(n, mean = mean1, sd = sd1) y <- rnorm(n, mean = mean2, sd = sd2) p.values[i] <- t.test(x, y)$p.value } hist(p.values, breaks=seq(0, 1, sig.level)) power <- sum(p.values <= sig.level) / trial.count return(power) }