# Here, we do a power analysis for a two-sample t test with
# true means of 100 and 110 and a standard deviation of 15.
# We assume that we have a total sample size of 25 (12 in one
# group and 15 in the other). 

# First, we calculate the noncentrality parameter:

> delta <- (110-100)/15 * sqrt((12*13/25))
> delta
[1] 1.665333

# Here's the critical value of t:

> tcrit <- qt(.975, 23)
> tcrit
[1] 2.068658

# So the power is:

> pt(-tcrit, 23, delta) + (1-pt(tcrit, 23, delta))
[1] 0.3581033

# (Not so great.)

# Next, we consider power analysis for an ANOVA in which
# five means are spread over a +- 1 sd range. (The details
# for the noncentrality parameter are in today's Powerpoint).

# We want to know what n we will need in each group in order
# to attain power of .90 (assuming an alpha level of .05).
# Here are some possible values of n:

> n <- 3:15

# For each n, the noncentrality parameter will differ:

> 2.5*n -> lambda
> cbind(n,lambda)
       n lambda
 [1,]  3    7.5
 [2,]  4   10.0
 [3,]  5   12.5
 [4,]  6   15.0
 [5,]  7   17.5
 [6,]  8   20.0
 [7,]  9   22.5
 [8,] 10   25.0
 [9,] 11   27.5
[10,] 12   30.0
[11,] 13   32.5
[12,] 14   35.0
[13,] 15   37.5

# There are 5 groups, so the degrees of freedom for the numerator
# of the F statistic will always be 4; but the df for the denominator
# will vary as n changes:

> dfn <- 4
> dfd <- 5*(n-1)

# And, because the df in the denominator changes, the critical value
# of F will also change:

> fcrit <- qf(.95, dfn, dfd)
> cbind(n,lambda,fcrit)
       n lambda    fcrit
 [1,]  3    7.5 3.478050
 [2,]  4   10.0 3.055568
 [3,]  5   12.5 2.866081
 [4,]  6   15.0 2.758710
 [5,]  7   17.5 2.689628
 [6,]  8   20.0 2.641465
 [7,]  9   22.5 2.605975
 [8,] 10   25.0 2.578739
 [9,] 11   27.5 2.557179
[10,] 12   30.0 2.539689
[11,] 13   32.5 2.525215
[12,] 14   35.0 2.513040
[13,] 15   37.5 2.502656

# Here, we calculate the power for each sample size:

> power <- 1-pf(fcrit, dfn, dfd, lambda)
> cbind(n, power)
       n     power
 [1,]  3 0.3827862
 [2,]  4 0.5645058
 [3,]  5 0.7110347
 [4,]  6 0.8177333
 [5,]  7 0.8898304
 [6,]  8 0.9357913
 [7,]  9 0.9637413
 [8,] 10 0.9800839
 [9,] 11 0.9893256
[10,] 12 0.9944029
[11,] 13 0.9971226
[12,] 14 0.9985470
[13,] 15 0.9992781

# So to achieve power of .9, we would need 8 participants
# in each group.

 
# The one situation in which a post hoc power analysis makes
# sense is when we are planning a replication. Let's suppose
# that we want to replicate the Eysenck experiment. We can use
# the results of our original experiment to drive the power analysis.
# Here are the means:

> attach(Eysenck)
> tapply(Score, Group, mean)
adjective  counting   imagery    intent   rhyming 
     11.0       7.0      13.4      12.0       6.9 

# We also need the error mean square:

> anova(lm(Score~Group))
Analysis of Variance Table

Response: Score
          Df Sum Sq Mean Sq F value    Pr(>F)    
Group      4 351.52  87.880  9.0848 1.815e-05 ***
Residuals 45 435.30   9.673                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# The square root of the error mean square will be our standard deviation:

> sqrt(9.673)
[1] 3.110145
 
# (We used those numbers for a G*power analysis.)


# Here's how to read in the data for Homework Five:

> Hwk5 <- read.csv("http://faculty.ucmerced.edu/jvevea/classes/105/data/Hwk5.csv")
> head(Hwk5)
    Group Score
1 Control    17
2 Control    12
3 Control    14
4 Control    17
5 Control    12
6 Control    11
>