quick example

We start with the entire dataframe df. Then we select only the rows where the height is taller than 100. Then we group by homeworld, and compute the mean birth year. What we get is a new data.frame, showing the mean birth years for each homeworld.

This is a common refrain:

Dataframe %>% filter %>% group_by %>% summarise

library(dplyr)

new_df <- df %>%
            filter(height > 100) %>%
            group_by(homeworld) %>%
            summarise(mean_birth_year = mean(birth_year,na.rm=TRUE))

filter

We can filter the data by properties of particular columns

# make a new dataframe that only include rows where the height is greater than 120
new_df <- df %>%
        filter(height > 120)
# multiple filters are seperated by a comma
new_df <- df %>%
        filter(height > 120,
               height < 180,
               birth_year > 20)
# more examples of differen logical operators
new_df <- df %>%
        filter(gender == "male") # == equals identity (same as)
new_df <- df %>%
        filter(gender != "male") # != not equal to
new_df <- df %>%
        filter(eye_color %in% c("blue","yellow") == TRUE) # looks for matches to blue and yellow
# <= less than or equal to
# >= greater than or equal to
new_df <- df %>%
        filter(height >= 120,
               height <= 180) 
new_df <- df %>%
        filter(height > 120 & height < 180) # &  AND
new_df <- df %>%
        filter(skin_color == "fair" | skin_color == "gold") # | OR

group_by and summarise

group_by() let’s us grab parts of the data based on the levels in the column

summarise() applies a function to each of the groups that are created

# counts the number of names, for each hair color
new_df <- df %>%
          group_by(hair_color) %>%
          summarise(counts=length(name))
# counts names, for each combination of hair and eye color
new_df <- df %>%
          group_by(hair_color,eye_color) %>%
          summarise(counts=length(name))

focus on summarise

summarise is very powerful. Using summarise we can apply any function we want to each of the groups. This includes intrinsic R functions, and functions of our own design. And, we can add as many as we want. What we get back are new dataframes with columns for each group, and new columns with variables containing the data we want from the analysis

new_df <- df %>%
          group_by(hair_color,eye_color) %>%
          summarise(mean_years = mean(birth_year,na.rm=TRUE),
                    sd_years = sd(birth_year,na.rm=TRUE),
                    counts = length(name))

Mutate

Use mutate to change or add a column

# change numbers in the height column by subtracting 100
library(dplyr)
new_df <- df %>%
            mutate(height=height-100)
# make a new column dividing height by mass
new_df <- df %>%
            mutate(hm_ratio = height/mass)

Select

Use select to select columns of interest and return a dataframe only with those columns

new_df <- df %>%
            select(name,height,mass)

Star wars questions and answers

1. use dplyr to find how many movies each character appeared in. Return a table that lists the names, and the number of films

# two ways to do the same thing
new_df <- df %>% 
            select(name,films) %>%
            group_by(name) %>%
            mutate(films = length(unlist(films)))
new_df <- df %>% 
            select(name,films) %>%
            group_by(name) %>%
            summarise(films = length(unlist(films)))

2. How many characters are in each movie?

• dplyr isn’t always necessary for every question

table(unlist(df$films))

## 
##              A New Hope    Attack of the Clones      Return of the Jedi 
##                      18                      40                      20 
##     Revenge of the Sith The Empire Strikes Back       The Force Awakens 
##                      34                      16                      11 
##      The Phantom Menace 
##                      34

3. What is the mean height and mass of characters with blue eyes?

new_df <- df %>%
            filter(eye_color=="blue") %>%
            summarise(mean_height = mean(height,na.rm=TRUE),
                      mean_mass = mean(mass,na.rm=TRUE))

Other helpful libraries for reading in files

1. libary readxl let’s you read in excel files
2. library googlesheets let’s you read in google spreadsheets
3. function scan is a powerful and all purpose text input function, often helpful in very messy situations where you want to read in line-by-line
4. load for R data files
5. there are several read. functions for specific situations.
6. library jsonlite for json data structures
7. and many more, this is really something that you will learn more about on a case-by-case basis as you have to deal with particular data formats.

one sample t-test

Set mu to test against a population value

x <- c(1,4,3,2,3,4,3,2,3,2,3,4)
#non-directional
t.test(x, mu = 2)

## 
##  One Sample t-test
## 
## data:  x
## t = 3.0794, df = 11, p-value = 0.01048
## alternative hypothesis: true mean is not equal to 2
## 95 percent confidence interval:
##  2.237714 3.428952
## sample estimates:
## mean of x 
##  2.833333

t.test(x, mu = 2, alternative = "two.sided") # the default

## 
##  One Sample t-test
## 
## data:  x
## t = 3.0794, df = 11, p-value = 0.01048
## alternative hypothesis: true mean is not equal to 2
## 95 percent confidence interval:
##  2.237714 3.428952
## sample estimates:
## mean of x 
##  2.833333

# directional
t.test(x, mu = 2, alternative = "greater")

## 
##  One Sample t-test
## 
## data:  x
## t = 3.0794, df = 11, p-value = 0.005241
## alternative hypothesis: true mean is greater than 2
## 95 percent confidence interval:
##  2.34734     Inf
## sample estimates:
## mean of x 
##  2.833333

t.test(x, mu = 2, alternative = "less")

## 
##  One Sample t-test
## 
## data:  x
## t = 3.0794, df = 11, p-value = 0.9948
## alternative hypothesis: true mean is less than 2
## 95 percent confidence interval:
##      -Inf 3.319326
## sample estimates:
## mean of x 
##  2.833333

t = mean/SEM SEM tells you how things vary

paired-sample t-test

Set paired=TRUE

x <- c(1,4,3,2,3,4,3,2,3,2,3,4)
y <- c(3,2,5,4,3,2,3,2,1,2,2,3)
#non-directional
t.test(x, y, paired=TRUE)

## 
##  Paired t-test
## 
## data:  x and y
## t = 0.37796, df = 11, p-value = 0.7126
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.8038766  1.1372099
## sample estimates:
## mean of the differences 
##               0.1666667

t.test(x, y, paired=TRUE, alternative = "two.sided") # the default

## 
##  Paired t-test
## 
## data:  x and y
## t = 0.37796, df = 11, p-value = 0.7126
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.8038766  1.1372099
## sample estimates:
## mean of the differences 
##               0.1666667

# directional
t.test(x, y, paired=TRUE, alternative = "greater")

## 
##  Paired t-test
## 
## data:  x and y
## t = 0.37796, df = 11, p-value = 0.3563
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -0.6252441        Inf
## sample estimates:
## mean of the differences 
##               0.1666667

t.test(x, y, paired=TRUE, alternative = "less")

## 
##  Paired t-test
## 
## data:  x and y
## t = 0.37796, df = 11, p-value = 0.6437
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##       -Inf 0.9585774
## sample estimates:
## mean of the differences 
##               0.1666667

Independent sample t-test

Note: Default is Welch test, set var.equal=TRUE for t-test

x <- c(1,4,3,2,3,4,3,2,3,2,3,4)
y <- c(3,2,5,4,3,2,3,2,1,2,2,3)
#non-directional
t.test(x, y, var.equal=TRUE)

## 
##  Two Sample t-test
## 
## data:  x and y
## t = 0.40519, df = 22, p-value = 0.6893
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.6863784  1.0197117
## sample estimates:
## mean of x mean of y 
##  2.833333  2.666667

t.test(x, y, var.equal=TRUE, alternative = "two.sided") # the default

## 
##  Two Sample t-test
## 
## data:  x and y
## t = 0.40519, df = 22, p-value = 0.6893
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.6863784  1.0197117
## sample estimates:
## mean of x mean of y 
##  2.833333  2.666667

# directional
t.test(x, y, var.equal=TRUE, alternative = "greater")

## 
##  Two Sample t-test
## 
## data:  x and y
## t = 0.40519, df = 22, p-value = 0.3446
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -0.5396454        Inf
## sample estimates:
## mean of x mean of y 
##  2.833333  2.666667

t.test(x, y, var.equal=TRUE, alternative = "less")

## 
##  Two Sample t-test
## 
## data:  x and y
## t = 0.40519, df = 22, p-value = 0.6554
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##       -Inf 0.8729787
## sample estimates:
## mean of x mean of y 
##  2.833333  2.666667

printing

library(broom)
t_results <- t.test(x, y, var.equal=TRUE)
knitr::kable(tidy(t_results))

writing

The contents of R variables can be written into R Markdown documents.

t(22) =0.41, p = 0.689

# write your own function
report_t <- function(x){
  return(paste(c("t(",round(x$parameter, digits=2),") ",
        "= ",round(x$statistic, digits=2),", ",
        "p = ",round(x$p.value, digits=3)), collapse=""))
}
report_t(t_results)

## [1] "t(22) = 0.41, p = 0.689"

t_results

## 
##  Two Sample t-test
## 
## data:  x and y
## t = 0.40519, df = 22, p-value = 0.6893
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.6863784  1.0197117
## sample estimates:
## mean of x mean of y 
##  2.833333  2.666667

The results of my t-test were t(22) = 0.41, p = 0.689.

between-subjects 1 factor

# run anova
aov_out <- aov(DV~Factor1, all_data)
# summary
sum_out <- summary(aov_out)
# kable xtable printing
library(xtable)
knitr::kable(xtable(summary(aov_out)))

# print means
print(model.tables(aov_out,"means"), format="markdown")

## Tables of means
## Grand mean
##            
## -0.1673187 
## 
##  Factor1 
## Factor1
##       A       B 
## -0.3231 -0.0115

sum_out[[1]]$`F value`[1]

## [1] 1.071932

#kable function good for making dataframes into pretty table

between-subjects 2 factor

# run anova
aov_out <- aov(DV~Factor1*Factor2, all_data)
# summary
summary(aov_out)

##                 Df Sum Sq Mean Sq F value Pr(>F)
## Factor1          1   0.97  0.9712   1.093  0.303
## Factor2          1   1.39  1.3881   1.562  0.219
## Factor1:Factor2  1   1.05  1.0498   1.181  0.284
## Residuals       36  31.99  0.8887

# kable xtable printing
library(xtable)
knitr::kable(xtable(summary(aov_out)))

# print means
print(model.tables(aov_out,"means"), format="markdown")

## Tables of means
## Grand mean
##            
## -0.1673187 
## 
##  Factor1 
## Factor1
##       A       B 
## -0.3231 -0.0115 
## 
##  Factor2 
## Factor2
##       1       2 
##  0.0190 -0.3536 
## 
##  Factor1:Factor2 
##        Factor2
## Factor1 1       2      
##       A  0.0251 -0.6714
##       B  0.0128 -0.0358

repeated-measures 1 factor

# run anova
aov_out <- aov(DV~Factor1 + Error(Subjects/Factor1), all_data)
# summary
summary(aov_out)

## 
## Error: Subjects
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  9  4.337  0.4819               
## 
## Error: Subjects:Factor1
##           Df Sum Sq Mean Sq F value Pr(>F)
## Factor1    1  0.971  0.9712   0.496  0.499
## Residuals  9 17.620  1.9578               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 20  12.47  0.6237

# kable xtable printing
library(xtable)
knitr::kable(xtable(summary(aov_out)))

# print means
print(model.tables(aov_out,"means"), format="markdown")

## Tables of means
## Grand mean
##            
## -0.1673187 
## 
##  Factor1 
## Factor1
##       A       B 
## -0.3231 -0.0115

repeated-measures 2 factor

# run anova
aov_out <- aov(DV~Factor1*Factor2 + Error(Subjects/(Factor1*Factor2)), all_data)
# summary
summary(aov_out)

## 
## Error: Subjects
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  9  4.337  0.4819               
## 
## Error: Subjects:Factor1
##           Df Sum Sq Mean Sq F value Pr(>F)
## Factor1    1  0.971  0.9712   0.496  0.499
## Residuals  9 17.620  1.9578               
## 
## Error: Subjects:Factor2
##           Df Sum Sq Mean Sq F value Pr(>F)
## Factor2    1  1.388  1.3881   1.726  0.221
## Residuals  9  7.239  0.8043               
## 
## Error: Subjects:Factor1:Factor2
##                 Df Sum Sq Mean Sq F value Pr(>F)  
## Factor1:Factor2  1  1.050  1.0498   3.378 0.0992 .
## Residuals        9  2.797  0.3108                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# kable xtable printing
library(xtable)
knitr::kable(xtable(summary(aov_out)))

# print means
print(model.tables(aov_out,"means"), format="markdown")

## Tables of means
## Grand mean
##            
## -0.1673187 
## 
##  Factor1 
## Factor1
##       A       B 
## -0.3231 -0.0115 
## 
##  Factor2 
## Factor2
##       1       2 
##  0.0190 -0.3536 
## 
##  Factor1:Factor2 
##        Factor2
## Factor1 1       2      
##       A  0.0251 -0.6714
##       B  0.0128 -0.0358

papaja

#library(papaja)
#apa_stuff <- apa_print.aov(aov_out)

#The main effect for factor 1 was, `r apa_stuff$statistic$Factor1`. The main effect for factor 2 was, `r apa_stuff$statistic$Factor2`. The interaction was, `r apa_stuff$statistic$Factor1_Factor2`

pairwise t-tests

Performs all possible t-tests between the levels of a grouping factor. Can implement different corrections for multiple comparisons.

df <- airquality #loads the airquality data from base R
df$Month <- as.factor(df$Month)
# default uses holm correction
pairwise.t.test(df$Ozone,df$Month)

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  df$Ozone and df$Month 
## 
##   5       6       7       8      
## 6 1.00000 -       -       -      
## 7 0.00026 0.05113 -       -      
## 8 0.00019 0.04987 1.00000 -      
## 9 1.00000 1.00000 0.00488 0.00388
## 
## P value adjustment method: holm

# use bonferroni correction
pairwise.t.test(df$Ozone,df$Month, p.adj="bonf")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  df$Ozone and df$Month 
## 
##   5       6       7       8      
## 6 1.00000 -       -       -      
## 7 0.00029 0.10225 -       -      
## 8 0.00019 0.08312 1.00000 -      
## 9 1.00000 1.00000 0.00697 0.00485
## 
## P value adjustment method: bonferroni

tukey HSD test

notes: the warpbreaks data.frame should be available in base R, you do not need to load it to be able to use it.

Doesn’t work with repeated measures ANOVAs.

summary(fm1 <- aov(breaks ~ wool + tension, data = warpbreaks))

##             Df Sum Sq Mean Sq F value  Pr(>F)   
## wool         1    451   450.7   3.339 0.07361 . 
## tension      2   2034  1017.1   7.537 0.00138 **
## Residuals   50   6748   135.0                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(fm1, "tension", ordered = TRUE)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
##     factor levels have been ordered
## 
## Fit: aov(formula = breaks ~ wool + tension, data = warpbreaks)
## 
## $tension
##          diff        lwr      upr     p adj
## M-H  4.722222 -4.6311985 14.07564 0.4474210
## L-H 14.722222  5.3688015 24.07564 0.0011218
## L-M 10.000000  0.6465793 19.35342 0.0336262

Fishers LSD

Explanation of Fisher’s least significant difference test. https://www.utd.edu/~herve/abdi-LSD2010-pretty.pdf

library(agricolae)
data(sweetpotato)
model <- aov(yield~virus, data=sweetpotato)
out   <- LSD.test(model,"virus", p.adj="bonferroni")

Null for a t-test

# samples A and B come from the same normal distribution
A <- rnorm(n=10,mean=10, sd=5)
B <- rnorm(n=10,mean=10, sd=5)
# the pvalue for this one pretend simulation
t.test(A,B,var.equal=TRUE)$p.value

## [1] 0.3711599

# running the simulation
# everytime we run this function we do one simulated experiment and return the p-value
sim_null <- function(){
  A <- rnorm(n=10,mean=10, sd=5)
  B <- rnorm(n=10,mean=10, sd=5)
  return(t.test(A,B,var.equal=TRUE)$p.value)
}
# use replicate to run the sim many times
outcomes <- replicate(1000,sim_null())
# plot the null-distribution of p-values
hist(outcomes)

# proportion of simulated experiments had a p-value less than .05
length(outcomes[outcomes<.05])/1000

## [1] 0.046

We ran the above simulation 1000 times. By definition, we should get approximately 5% of the simulations returning a p-value less than .05. If we increase the number of simulations, then we will get a more accurate answer that converges on 5% every time.

pwr package

Here is an example of using the pwr package to find the power for an independent sample t-test, with n=10, to detect an effect-size of 1. This answer is similar to our simulations answer. The simulation would converge on this answer if we increased the number of simulations.

library(pwr)
pwr.t.test(n=10,
           d=1,
           sig.level=.05,
           type="two.sample",
           alternative="two.sided")

## 
##      Two-sample t test power calculation 
## 
##               n = 10
##               d = 1
##       sig.level = 0.05
##           power = 0.5620066
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

As I mentioned there are many functions for directly computing power in R. Feel free to use them. In this class, we learn how to use simulation to conduct power analyses. This can be a redundant approach that is not necessary, given there are other functions we can use. Additionally, we will not get exact solutions (but approximate ones). Nevertheless, the existing power functions can be limited and may not apply to your design. The simulation approach can be extended to any design. Learning how to run the simulations will also improve your statistical sensibilities, and power calculations will become less of a black box.

Two more things before moving onto simulation: power-curves, and sample-size planning.

power-curves

A design’s power is in relation to the true effect-size. The same design has different levels of power to detect different sized effects. Let’s make a power curve to see the power of a t-test for independent samples (n=10) to detect a range of effect-sizes

effect_sizes <- seq(.1,2,.1)
power <- sapply(effect_sizes, 
          FUN = function(x) {
            pwr.t.test(n=10,
            d=x,
            sig.level=.05,
            type="two.sample",
            alternative="two.sided")$power})
plot_df <- data.frame(effect_sizes,power)
library(ggplot2)
ggplot(plot_df, aes(x=effect_sizes,
                    y=power))+
  geom_point()+
  geom_line()

This power curve applies to all independent-sample t-tests with n=10. It is a property, or fact about those designs. Every design has it’s own power curve. The power curve shows us what should happen (on average), when the true state of the world involves effects of different sizes.

If you do not know the power-curve for your design, then you do not know how sensitive your design is for detecting effects of particular sizes. You might accidentally be using an under powered design, with only a very small chance of detecting an effect of a size you are interested in.

If you do know the power-curve for your design, you will be in a better position to plan your experiment, for example by modifying the number of subjects that you run.

Sample-size planning

Here is one way to plan for the number of subjects that you need to find an effect of interest.

1. Establish a smallest effect-size of interest
2. Create a curve showing the power of your design as a function of number of subjects to detect the smallest effect-size of interest.

It’s not clear how you establish a smallest effect-size of interest. But let’s say you are interested in detecting an effect of at least d = .2. This means that two conditions would differ by at least a .2 standard deviation shift. If you find something smaller than that, let’s say you wouldn’t care about it because it wouldn’t be big enough for you to care. How many subjects do you need to have a high powered design, one that would almost always reject the null-hypothesis?

This is for an independent samples t-test:

num_subjects <- seq(10,1000,10)
power <- sapply(num_subjects, 
          FUN = function(x) {
            pwr.t.test(n=x,
            d=.2,
            sig.level=.05,
            type="two.sample",
            alternative="two.sided")$power})
plot_df <- data.frame(num_subjects,power)
library(ggplot2)
ggplot(plot_df, aes(x=num_subjects,
                    y=power))+
  geom_point()+
  geom_line()

Well, it looks like you need many subjects to have high power. For example, if you want to detect the effect 95% of the time, you would need around 650 subjects. It’s worth doing this kind of analysis to see if your design checks out. You don’t want to waste your time running an experiment that is designed to fail (even when the true effect is real).

Simulated power for a t-test

A power curve for n=10.

# function to run a simulated t-test
sim_power <- function(x){
  A <- rnorm(n=10,mean=0, sd=1)
  B <- rnorm(n=10,mean=(0+x), sd=1)
  return(t.test(A,B,var.equal=TRUE)$p.value)
}
# vector of effect sizes
effect_sizes <- seq(.1,2,.1)
# run simulation for each effect size 1000 times
power <- sapply(effect_sizes, 
          FUN = function(x) {
            sims <- replicate(1000,sim_power(x))
            sim_power <- length(sims[sims<.05])/length(sims)
            return(sim_power)})
# combine into dataframe
plot_df <- data.frame(effect_sizes,power)
# plot the power curve
ggplot(plot_df, aes(x=effect_sizes,
                    y=power))+
  geom_point()+
  geom_line()

In this case, there is no obvious benefit to computing the power-curve by simulation. The answer we get is similar to the answer we got before using the pwr package, but our simulation answer is more noisy. Why bother the simulation?

One answer to the why bother question is that you can simulate deeper aspects of the design and get more refined answers without having to work through the math.

Simulating cell-size

Many experimental designs involve multiple measurements, or trials, for each subject in each condition. How many trials should we require for each subject in each condition? Traditional power analysis doesn’t make it easy to answer this question. However, the power of a design will depend not only on the number subjects, but also the number trials used to estimate the mean for each subject in each condition.

Consider a simple Stroop experiment. The researcher is interested in measuring a Stroop effect of at least d=.1 (e.g., the difference between mean congruent trials is .1 standard deviations smaller than mean incongruent trials). How many subjects are required? And, how many trials should each subject perform in the congruent and incongruent conditions? Let’s use simulation to find out.

# function to run a simulated t-test
# nsubs sets number of subjects
# ntrials to change number of trials
# d sets effect size
# this is a paired sample test to model Stroop
sim_power <- function(nsubs,ntrials,d){
  A <- replicate(nsubs,mean(rnorm(n=ntrials,mean=0, sd=1)))
  B <- replicate(nsubs,mean(rnorm(n=ntrials,mean=d, sd=1)))
  return(t.test(A,B,paired=TRUE)$p.value)
}
# vectors for number of subjects and trials
n_subs_vector <- c(10,20,30,50)
n_trials_vector <- c(10,20,30,50,100)
# a loop to run all simulations
power <- c()
subjects <- c()
trials <- c()
i <- 0 # use this as a counter for indexing
for(s in n_subs_vector){
  for(t in n_trials_vector){
    i <- i+1
    sims <- replicate(1000,sim_power(s,t,.1))
    power[i] <- length(sims[sims<.05])/length(sims)
    subjects[i] <- s
    trials[i] <- t
  }
}
# combine into dataframe
plot_df <- data.frame(power,subjects,trials)
# plot the power curve
ggplot(plot_df, aes(x=subjects,
                    y=power,
                    group=trials,
                    color=trials))+
  geom_point()+
  geom_line()

# a vectorized version of the loop
# run simulation for each effect size 1000 times
# power <- outer(n_subs_vector,
#                n_trials_vector,
#                FUN = Vectorize(function(x,y) {
#                   sims <- replicate(100,sim_power(x,y))  
#                   sim_power <- length(sims[sims<.05])/length(sims)
#                   return(sim_power)
#                }))

To my eye, it looks like 30 subjects with 100 trials in each condition would give you a very high power to find a Stroop effect of d=.1.

One-Way ANOVA

Let’s extend our simulation based approach to the one-way ANOVA. Let’s assume a between-subjects design, with one factor that has four levels: A, B, C, and D. There will be 20 subjects in each group. What is the power curve for this design to detect effects of various size? Immediately, the situation becomes complicated, there are numerous ways that the means for A, B, C, and D could vary. Let’s assume the simplest case, three of them are the same, and one of them is different by some amount of standard deviations. We will compute the main effect, and report the proportion of significant experiments as we increase the effect size for the fourth group.

# function to run a simulated t-test
sim_power_anova <- function(x){
  A <- rnorm(n=20,mean=0, sd=1)
  B <- rnorm(n=20,mean=0, sd=1)
  C <- rnorm(n=20,mean=0, sd=1)
  D <- rnorm(n=20,mean=(0+x), sd=1)
  df <- data.frame(condition = as.factor(rep(c("A","B","C","D"),each=20)),
                   DV = c(A,B,C,D))
  aov_results <- summary(aov(DV~condition,df))
  #return the pvalue
  return(aov_results[[1]]$`Pr(>F)`[1])
}
# vector of effect sizes
effect_sizes <- seq(.1,2,.1)
# run simulation for each effect size 1000 times
power <- sapply(effect_sizes, 
          FUN = function(x) {
            sims <- replicate(1000,sim_power_anova(x))
            sim_power <- length(sims[sims<.05])/length(sims)
            return(sim_power)})
# combine into dataframe
plot_df <- data.frame(effect_sizes,power)
# plot the power curve
ggplot(plot_df, aes(x=effect_sizes,
                    y=power))+
  geom_point()+
  geom_line()

It looks like this design (1 factor, between-subjects, 20 subjects per group) has high power to detect an effect of d=1, specifically when one of the groups differs from the others by d=1.

However, most effects in psychology are smalls, d=.2 is very common. How, many subjects does this design require to have high power (let’s say above .95, although most people use .8) to detect that small effect?

sim_power_anova <- function(x){
  A <- rnorm(n=x,mean=0, sd=1)
  B <- rnorm(n=x,mean=0, sd=1)
  C <- rnorm(n=x,mean=0, sd=1)
  D <- rnorm(n=x,mean=.2, sd=1)
  df <- data.frame(condition = as.factor(rep(c("A","B","C","D"),each=x)),
                   DV = c(A,B,C,D))
  aov_results <- summary(aov(DV~condition,df))
  #return the pvalue
  return(aov_results[[1]]$`Pr(>F)`[1])
}
# vector of effect sizes
subjects <- seq(10,1000,50)
# run simulation for each effect size 1000 times
power <- sapply(subjects, 
          FUN = function(x) {
            sims <- replicate(1000,sim_power_anova(x))
            sim_power <- length(sims[sims<.05])/length(sims)
            return(sim_power)})
# combine into dataframe
plot_df <- data.frame(subjects,power)
# plot the power curve
ggplot(plot_df, aes(x=subjects,
                    y=power))+
  geom_point()+
  geom_line()

The simulations suggests we need about 560 subjects in each group to have power .95 to detect the effect (d=.2). That’s a total of 2240 subjects. Reality can be surprising when it comes to power analysis. It is better to be surprised about your design before you run your experiment, not after.

Correlation between traits and behavior

A common research strategy is to measure putative correlations between people’s traits and their behavior. For example, a researcher might prepare a questionnaire with several questions. Perhaps the questions are about their political views, or about their life satisfaction, or anything else. The research might ask some group of people to answer the questions, and to also perform some task. At the end of the experiment, they might ask if different kinds of people (measured by the questionnaire) perform differently on the behavioral measure. For example, we might have people answer several questions about their openness to new experiences, we might also measure their performance on a working memory task. A research question might be, do people who are more open to experience perform better on a working memory task? The answer would lie in a correlation between the answers on the questions, and the performance on the task. Let’s simulate this kind of situation, say for 20 subjects. Each subject answers 20 questions, and they perform a behavioral task. At the end of the experiment, we correlate the answers on each question, with the performance on the task. What kind of correlations would we expect by chance alone? Let’s say the participants are all random robots. They answer each question randomly, and their performance on the behavioral task is sampled randomly from some normal distribution. If everything is random, shouldn’t we expect to find no correlations?

Details:

Let’s assume that each question involves a likert scale from 1 to 7. Each person randomly picks a number from 1 to 7 for each question. Let’s assume performance on the behavioral task is sampled from a normal distribution with mean = 0, and sd = 1.

# get 20 random answers for all 20 subjects and 20 questions
# columns will be individual subjects 1 to 20
# rows will be questions 1 to 20
questionnaire <- matrix(sample(1:7,20*20, replace=TRUE),ncol=20)
# get 20 measures of performance on behavioral task
behavioral_task <- rnorm(20,0,1)
# correlate behavior with each question
save_correlations <-c()
for(i in 1:20){
  save_correlations[i] <- cor(behavioral_task,questionnaire[i,])
}
# show histogram of 20 correlations
hist(save_correlations)

The histogram shows that a range of correlations between individual questions and behavior can emerge just by chance alone. If you run the above code a few times, you will see that the histogram changes a bit because of random chance.

Oftentimes researchers might not know which question on the questionnaire is the best question. That is, the one that best correlates with behavior. Consider a researcher who computes all of the correlations between each question and behavior, and then chooses the question with highest correlation (positive or negative) as the best question. After all, it has the highest correlation. After choosing this question, they might suggest that behavior is strongly correlated with how people answer this question.

Let’s try to find out by simulation what kinds of large correlations can occur just by chance alone. We will run the above many times, and each time we will save the absolute value of the largest correlation between a question and behavior. Just how large can these correlations be by chance alone?

save_max <- c()
for( j in 1:10000){
  questionnaire <- matrix(sample(1:7,20*20, replace=TRUE),ncol=20)
  behavioral_task <- rnorm(20,0,1)
  
  save_correlations <-c()
  for(i in 1:20){
    save_correlations[i] <- cor(behavioral_task,questionnaire[i,])
  }
  save_max[j] <- max(abs(save_correlations))
}
hist(save_max)

The simulation shows that chance alone in this situation can produce very large correlations, as large as .8 or .9 (although not very often).

The situation changes somewhat if many more subjects are run. Let’s do the same as above, but run 200 subjects, rather than 20.

save_max <- c()
for( j in 1:10000){
  questionnaire <- matrix(sample(1:7,200*20, replace=TRUE),ncol=200)
  behavioral_task <- rnorm(200,0,1)
  
  save_correlations <-c()
  for(i in 1:20){
    save_correlations[i] <- cor(behavioral_task,questionnaire[i,])
  }
  save_max[j] <- max(abs(save_correlations))
}
hist(save_max)

Now, chance doesn’t do much better than .25.