# Statistical power calculation in R using the package "pwr":
# pwr: https://cran.r-project.org/web/packages/pwr/index.html 

# If pwr is not installed yet: 
# install.packages("pwr") 

library(pwr)

# For t-test power calculation (equal sample sizes): 
# pwr.t.test(n = , d = , sig.level = , power = , type = c("two.sample", "one.sample", "paired"))
# d is for effect size 
# Cohen suggests that d values of 0.2, 0.5, and 0.8 represent small, medium, and large effect sizes, respectively.  
# For Cohen's d value: https://www.rdocumentation.org/packages/DescTools/versions/0.99.19/topics/CohenD 
pwr.t.test(n = 30, d = 0.5, sig.level = 0.05, type = "two.sample")

# For t-test power calculation (unequal sample sizes): 
# pwr.t2n.test(n1 = , n2= , d = , sig.level =, power = )
pwr.t2n.test(n1 = 30, n2= 50, d = 0.5, sig.level = 0.05)

# For a one-way analysis of variance power calculation: 
# pwr.anova.test(k = , n = , f = , sig.level = , power = )
# k is the number of groups, n is the common sample size in each group, and f is effect size
# Cohen suggests that f values of 0.1, 0.25, and 0.4 represent small, medium, and large effect sizes respectively
pwr.anova.test(k = 4, n = 12, f = 0.25, sig.level = 0.05)

# For correlation coefficients power analysis: 
# pwr.r.test(n = , r = , sig.level = , power = )
# n is the sample size and r is the correlation coefficient
# Cohen suggests that r values of 0.1, 0.3, and 0.5 represent small, medium, and large effect sizes, respectively
pwr.r.test(n = 50, r = 0.3, sig.level = 0.05) 

# When comparing two proportions from equal size samples, to calculate the power of the test:
# pwr.2p.test(h = , n = , sig.level = , power = )
# h is the effect size and n is the common sample size in each group
# Cohen suggests that h values of 0.2, 0.5, and 0.8 represent small, medium, and large effect sizes, respectively
pwr.2p.test(h = 0.5, n = 50, sig.level = 0.05)

# For unequal sample sizes:
# pwr.2p2n.test(h = , n1 = , n2 = , sig.level = , power = ) 
pwr.2p2n.test(h = 0.5, n1 = 40, n2 = 60, sig.level = 0.05)   # Example

# For chi-square tests:
# pwr.chisq.test(w = , N = , df = , sig.level = , power = )
# w is the effect size, N is the total sample size, and df is the degrees of freedom
# Cohen suggests that w values of 0.1, 0.3, and 0.5 represent small, medium, and large effect sizes, respectively
pwr.chisq.test(w = 0.3, N = 100, df = 1, sig.level = 0.05)   # Example
# An alternative method is provided in DescTools::power.chisq.test(). SEE: https://www.rdocumentation.org/packages/DescTools/versions/0.99.19/topics/power.chisq.test 

# EXAMPLES: Power calculations from Quick-R:
# For a one-way ANOVA comparing 5 groups, calculate the
# sample size needed in each group to obtain a power of
# 0.80, when the effect size is moderate (0.25) and a
# significance level of 0.05 is employed.
pwr.anova.test(k=5,f=.25,sig.level=.05,power=.8)

# What is the power of a one-tailed t-test, with a
# significance level of 0.01, 25 people in each group,
# and an effect size equal to 0.75?
pwr.t.test(n=25,d=0.75,sig.level=.01,alternative="greater")

# Using a two-tailed test proportions, and assuming a
# significance level of 0.01 and a common sample size of
# 30 for each proportion, what effect size can be detected
# with a power of .75?
pwr.2p.test(n=30,sig.level=0.01,power=0.75)

# Power calculation for GLM using the R package pwr
# https://www.statmethods.net/stats/power.html  
# https://data-se.netlify.com/2018/07/24/power-calculation-for-the-general-linear-model 

# For more examples, see: https://cran.r-project.org/web/packages/pwr/vignettes/pwr-vignette.html 

# Creating Power or Sample Size Plots
# Plot sample size curves for detecting correlations of various sizes

# range of correlations
r <- seq(.1, .5, .01)
nr <- length(r)

# power values
p <- seq(.4, .9, .1)
np <- length(p)

# obtain sample sizes
samsize <- array(numeric(nr*np), dim=c(nr,np))
for (i in 1:np){
  for (j in 1:nr){
    result <- pwr.r.test(n = NULL, r = r[j],
    sig.level = .05, power = p[i],
    alternative = "two.sided")
    samsize[j,i] <- ceiling(result$n)
  }
}

# set up graph
xrange <- range(r)
yrange <- round(range(samsize))
colors <- rainbow(length(p))
plot(xrange, yrange, type="n",
  xlab="Correlation Coefficient (r)",
  ylab="Sample Size (n)" )

# add power curves
for (i in 1:np){
  lines(r, samsize[,i], type="l", lwd=2, col=colors[i])
}

# add annotation (grid lines, title, legend)
abline(v=0, h=seq(0,yrange[2],50), lty=2, col="grey89")
abline(h=0, v=seq(xrange[1],xrange[2],.02), lty=2,
   col="grey89")
title("Sample Size Estimation for Correlation Studies\n
  Sig=0.05 (Two-tailed)")
legend("topright", title="Power", as.character(p),
   fill=colors)


# t-test power curve plot:
# FROM: https://www.r-exercises.com/2017/11/11/power-analysis-tutorial

sample_sizes = seq(5, 300, by=5)
power_t_test = pwr.t.test(n = sample_sizes , d = 0.3 , sig.level = 0.05, power = , type = "two.sample")

par(mfrow=c(1,1),oma=c(0,0,0,0),mar=c(4.5,4.5,0.5,0.5))
plot(sample_sizes,power_t_test$power, type='b', pch=16, las=1,
     ylab='statistical power', xlab='sample size',
     cex.lab=1.4, cex.axis=1.2, ylim = c(0,1), col = "blue")
abline(h = 0.80, lty=2, lwd=2,col='red')


# For effect sizes, see: https://cran.r-project.org/web/packages/pwr/vignettes/pwr-vignette.html

# To find out an effect size (ES):
cohen.ES(test = "r", size = "medium")
cohen.ES(test = "t", size = "medium")
cohen.ES(test = "anov", size = "medium")
cohen.ES(test = "p", size = "medium")
cohen.ES(test = "chisq", size = "medium")
cohen.ES(test = "f2", size = "medium")
# 'arg(test)' should be one of “p”, “t”, “r”, “anov”, “chisq”, “f2”
# 'arg(size)' should be one of “small”, “medium”, “large”


# List of functions in pwr:
# function 	power calculations for
# pwr.2p.test 	two proportions (equal n)
# pwr.2p2n.test 	two proportions (unequal n)
# pwr.anova.test 	balanced one way ANOVA
# pwr.chisq.test 	chi-square test
# pwr.f2.test 	general linear model
# pwr.p.test 	proportion (one sample)
# pwr.r.test 	correlation
# pwr.t.test 	t-tests (one sample, 2 sample, paired)
# pwr.t2n.test 	t-test (two samples with unequal n) 
# FROM: https://www.statmethods.net/stats/power.html