Power analysis (and other stuff)!

# REFRESHER: POWER is the probability of correctly rejecting a null hypothesis that is actually false.
# P(rejecting NULL hypothesis|NULL hypothesis is false)
# In other words, its the prob. of finding an effect that is really present.
# Power = 1 - P(Type II error)

# Power, effect size, sample size, and the significance level are inter-related, and if you know 3 of these
# quantities you can calculate the fourth (exciting eh?). Let's get to it.

#First, we'll need to install the "pwr" package.
install.packages("pwr")
library(pwr)

# Calculating Cohen's effect size (Rules of Thumb) in case we werent able to remember.
# select the stats test of interest and one effect size from:
cohen.ES(test = c("p", "t", "r", "anov", "chisq", "f2"), size = c("small", "medium", "large"))
# to determine effect size according to Cohen's Rules of Thumb.
# For example, a large effect for a multiple regression would be:
cohen.ES(test = c("f2"), size = c("large"))
# whereas a small effect for anova would be:
cohen.ES(test = c("chisq"), size = c("medium"))

#POWER ANALYSIS FUNCTIONS!

# T-TESTS!
# Function:
pwr.t.test(n = , d = , sig.level = , power = , type = c("two.sample", "one.sample", "paired"),
alternative= c("two.sided", "less", "greater"))
# Remember, power, effect size, sample size, and the significance level are related so we can play around
# with the information we have to determine the missing quantity.

#For example, what is the power for a two-sided ind. samples t-test given n=25 and a medium effect size?
pwr.t.test(n=25, d=.5, sig.level= .05)

#Another example: What is the sample size needed to achieve a pwr=.80 and a small effect size?
pwr.t.test(d=.2, sig.level= .05, power=.8)
# Note: if you add "$n" to the end of the command only the necessary sample size will be displayed.
pwr.t.test(d=.2, sig.level= .05, power=.8)$n

# R can also calculate pwr for two samples with unequal sample sizes:
pwr.t2n.test(n1= 35, n2= 17, d=.8, sig.level= .05, alternative = "greater")

# ANOVA!
# what's Cohen's rules of thumb for anova effect size (i.e., "f")?
cohen.ES(test = c("anov"), size = c("small"))
cohen.ES(test = c("anov"), size = c("medium"))
cohen.ES(test = c("anov"), size = c("large"))

# One-way ANOVA pwr analysis function:
pwr.anova.test(k = , n = , f = , sig.level = , power = )
# where k is the number of groups and n is the number of obs. IN EACH GROUP.
pwr.anova.test(k =4, n=32, f=.3, sig.level =.05)

# More complex (e.g., factorial ANOVA, ANCOVA) pwr analysis?? --> to be continued.

# PROPORTIONS!
# Function:
pwr.2p.test(h = , n = , sig.level =, power = )              # 2 proportions with equal sample sizes
pwr.2p2n.test(h = , n1 = , n2 = , sig.level = , power = )   # 2 proportions with unequal sample sizes

# CORRELATION!
#Function:
pwr.r.test(n = , r = , sig.level = , power = )

# CHI-SQUARE
# Function:
pwr.chisq.test(w =, N = , df = , sig.level =, power = )
# where "w"= effect size and "n"=number of observations.
# Quick example:
pwr.chisq.test(w=.3, N=120, df=9, sig.level=.05)

# REGRESSION!
# Function:
pwr.f2.test(u =, v = , f2 = , sig.level = , power = )
# where "u"= numerator degrees of freedom (number of continuous variables + number of dummy codes - 1)
# and "v"=denominator (error) degrees of freedom

# For example, find the power for a multiple regression test with 2 continuous predictors and 1 categorical
# predictor (i.e. Marital status with k=3 so 3-1=2 dummy codes) that has a large effect size and a sample size of 30.
pwr.f2.test(u =3, v =30, f2 =.35, sig.level =.05)

# Generating a table of sample sizes
# What is the required sample size for an independent samples t-test with pwr=.80?
seq <- seq(.1,.8, .1)
FindN <- rep(0,8)
for(i in 1:8) FindN[i]=pwr.t.test(d=seq[i],power=.8,sig.level=.05,type="two.sample")$n
data.frame(d=seq , N=ceiling(FindN))
# can swap in other functions, just make sure arguments are relevant to test
seq <- seq(.1,.8, .1)
FindN <- rep(0,8)
for(i in 1:8) FindN[i]=pwr.anova.test(k=4, f=seq[i], power=.8,sig.level=.05)$n
data.frame(f=seq , N=ceiling(FindN))
# NOTE: "ceiling" is an argument for rounding numbers so the values are slightly different than if we
# computed the sample size based on specified values. For example...
pwr.anova.test(k =4,f=.3, sig.level =.05, power=.8)

# Finally, plot of sample size curves for detecting effect at various sample sizes.
# range of correlations
d <- seq(.1,.5,.01)
nd <- length(d)

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

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

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

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

# add title and legend
title("Sample Size Estimation for t-test")
legend("topright", title="Power", as.character(p),
fill=colors)