14-AppendixBrysbaert.Rmd

# Appendix 3: Comparison to Brysbaert {-}

## One-Way ANOVA

Now we will simply replicate the simulations of @brysbaert2019many, and compare those results to `Superpower`. Simulations to estimate the power of an ANOVA with three unrelated groups the effect between the two extreme groups is set to d = .4,  the effect for the third group is d = .4 (see below for other situations) we use the built-in aov-test command give sample sizes (all samples sizes are equal).

```{r brysbaert_loop_1, eval=FALSE}
# Simulations to estimate the power of an ANOVA 
#with three unrelated groups
# the effect between the two extreme groups is set to d = .4, 
# the effect for the third group is d = .4 
#(see below for other situations)
# we use the built-in aov-test command
# give sample sizes (all samples sizes are equal)
N = 90
# give effect size d
d1 = .4 # difference between the extremes
d2 = .4 # third condition goes with the highest extreme
# give number of simulations
nSim = nsims
# give alpha levels
# alpha level for the omnibus ANOVA
alpha1 = .05 
#alpha level for three post hoc one-tail t-tests Bonferroni correction
alpha2 = .05 
```

```{r eval=FALSE}
# create vectors to store p-values
p1 <- numeric(nSim) #p-value omnibus ANOVA
p2 <- numeric(nSim) #p-value first post hoc test
p3 <- numeric(nSim) #p-value second post hoc test
p4 <- numeric(nSim) #p-value third post hoc test
pes1 <- numeric(nSim) #partial eta-squared
pes2 <- numeric(nSim) #partial eta-squared two extreme conditions
```

```{r eval=FALSE}
for (i in 1:nSim) {

x <- rnorm(n = N, mean = 0, sd = 1)
y <- rnorm(n = N, mean = d1, sd = 1)
z <- rnorm(n = N, mean = d2, sd = 1)
data = c(x, y, z)
groups = factor(rep(letters[24:26], each = N))
test <- aov(data ~ groups)
pes1[i] <- etaSquared(test)[1, 2]
p1[i] <- summary(test)[[1]][["Pr(>F)"]][[1]]
p2[i] <- t.test(x, y)$p.value
p3[i] <- t.test(x, z)$p.value
p4[i] <- t.test(y, z)$p.value
data = c(x, y)
groups = factor(rep(letters[24:25], each = N))
test <- aov(data ~ groups)
pes2[i] <- etaSquared(test)[1, 2]
}
```

```{r eval=FALSE}
# results are as predicted when omnibus ANOVA is significant,
# t-tests are significant between x and y plus x and z; 
# not significant between y and z
# printing all unique tests (adjusted code by DL)
sum(p1 < alpha1) / nSim
sum(p2 < alpha2) / nSim
sum(p3 < alpha2) / nSim
sum(p4 < alpha2) / nSim
mean(pes1)
mean(pes2)
```
```{r include=FALSE}
#p1
p1_1 
#p2
p2_1 
#p3
p3_1 
#4
p4_1 
#pes1
pes1_1 
#pes2
pes2_1 
```

### Three conditions replication

```{r}
K <- 3
mu <- c(0, 0.4, 0.4)
n <- 90
sd <- 1
r <- 0
design = paste(K, "b", sep = "")
```

```{r, message=FALSE, warning=FALSE}
design_result <- ANOVA_design(
  design = design,
  n = n,
  mu = mu,
  sd = sd,
  labelnames = c("factor1", "level1", "level2", "level3")
  )
```
```{r eval=FALSE, echo=TRUE}
simulation_result <- ANOVA_power(design_result, 
                                 alpha_level = alpha_level, 
                                 nsims = nsims,
                                 verbose = FALSE)

```

```{r echo=FALSE}
knitr::kable(simulation_result_2.11$main_results,
             caption = "Simulated ANOVA Result")%>%
  kable_styling(latex_options = "hold_position")
```


```{r}
exact_result <- ANOVA_exact(design_result,
                            alpha_level = alpha_level,
                            verbose = FALSE)
```


```{r echo=FALSE}
knitr::kable(exact_result$main_results,
             caption = "Exact ANOVA Result")%>%
  kable_styling(latex_options = "hold_position")
```


### Variation 1

```{r eval=FALSE}
# give sample sizes (all samples sizes are equal)
N = 145
# give effect size d
d1 = .4 #difference between the extremes
d2 = .0 #third condition goes with the highest extreme
# give number of simulations
nSim = nsims
# give alpha levels
#alpha level for the omnibus ANOVA
alpha1 = .05 
#alpha level for three post hoc one-tail t-test Bonferroni correction
alpha2 = .05 
```

```{r eval=FALSE}
# create vectors to store p-values
p1 <- numeric(nSim) #p-value omnibus ANOVA
p2 <- numeric(nSim) #p-value first post hoc test
p3 <- numeric(nSim) #p-value second post hoc test
p4 <- numeric(nSim) #p-value third post hoc test
pes1 <- numeric(nSim) #partial eta-squared
pes2 <- numeric(nSim) #partial eta-squared two extreme conditions
```

```{r eval=FALSE}
for (i in 1:nSim) {


x <- rnorm(n = N, mean = 0, sd = 1)
y <- rnorm(n = N, mean = d1, sd = 1)
z <- rnorm(n = N, mean = d2, sd = 1)
data = c(x, y, z)
groups = factor(rep(letters[24:26], each = N))
test <- aov(data ~ groups)
pes1[i] <- etaSquared(test)[1, 2]
p1[i] <- summary(test)[[1]][["Pr(>F)"]][[1]]
p2[i] <- t.test(x, y)$p.value
p3[i] <- t.test(x, z)$p.value
p4[i] <- t.test(y, z)$p.value
data = c(x, y)
groups = factor(rep(letters[24:25], each = N))
test <- aov(data ~ groups)
pes2[i] <- etaSquared(test)[1, 2]
}
```

```{r eval=FALSE}
# results are as predicted when omnibus ANOVA is significant, 
# t-tests are significant between x and y plus x and z; 
# not significant between y and z
# printing all unique tests (adjusted code by DL)
sum(p1 < alpha1) / nSim
sum(p2 < alpha2) / nSim
sum(p3 < alpha2) / nSim
sum(p4 < alpha2) / nSim
mean(pes1)
mean(pes2)
```
```{r include=FALSE}
#p1
p1_2 
#p2
p2_2 
#p3
p3_2 
#4
p4_2 
#pes1
pes1_2 
#pes2
pes2_2 
```

### Three conditions replication


```{r}
K <- 3
mu <- c(0, 0.4, 0.0)
n <- 145
sd <- 1
r <- 0
design = paste(K, "b", sep = "")
```

```{r, message=FALSE, warning=FALSE}
design_result <- ANOVA_design(
  design = design,
  n = n,
  mu = mu,
  sd = sd,
  labelnames = c("factor1", "level1", "level2", "level3")
  )
```
```{r eval=FALSE, echo=TRUE}
simulation_result <- ANOVA_power(design_result, 
                                 alpha_level = alpha_level, 
                                 nsims = nsims,
                                 verbose = FALSE)

```

```{r echo = FALSE}
knitr::kable(simulation_result_2.13$main_results,
             caption = "Simulated ANOVA Result")%>%
  kable_styling(latex_options = "hold_position")
```


```{r}
exact_result <- ANOVA_exact(design_result,
                            alpha_level = alpha_level,
                            verbose = FALSE)
```


```{r echo = FALSE}
knitr::kable(exact_result$main_results,
             caption = "Exact ANOVA Result")%>%
  kable_styling(latex_options = "hold_position")
```

### Variation 2


```{r eval=FALSE}
# give sample sizes (all samples sizes are equal)
N = 82
# give effect size d
d1 = .4 #difference between the extremes
d2 = .2 #third condition goes with the highest extreme
# give number of simulations
nSim = nsims
# give alpha levels
#alpha level for the omnibus ANOVA
alpha1 = .05 
#alpha level for three post hoc one-tail t-test Bonferroni correction
alpha2 = .05 
```

```{r eval=FALSE}
# create vectors to store p-values
p1 <- numeric(nSim) #p-value omnibus ANOVA
p2 <- numeric(nSim) #p-value first post hoc test
p3 <- numeric(nSim) #p-value second post hoc test
p4 <- numeric(nSim) #p-value third post hoc test
pes1 <- numeric(nSim) #partial eta-squared
```

```{r eval=FALSE}
for (i in 1:nSim) {
#for each simulated experiment

x <- rnorm(n = N, mean = 0, sd = 1)
y <- rnorm(n = N, mean = d1, sd = 1)
z <- rnorm(n = N, mean = d2, sd = 1)
data = c(x, y, z)
groups = factor(rep(letters[24:26], each = N))
test <- aov(data ~ groups)
pes1[i] <- etaSquared(test)[1, 2]
p1[i] <- summary(test)[[1]][["Pr(>F)"]][[1]]
p2[i] <- t.test(x, y)$p.value
p3[i] <- t.test(x, z)$p.value
p4[i] <- t.test(y, z)$p.value
data = c(x, y)
groups = factor(rep(letters[24:25], each = N))
test <- aov(data ~ groups)
pes2[i] <- etaSquared(test)[1, 2]
}
```

```{r eval=FALSE}
sum(p1 < alpha1) / nSim
sum(p2 < alpha2) / nSim
sum(p3 < alpha2) / nSim
sum(p4 < alpha2) / nSim
mean(pes1)
mean(pes2)
```
```{r include=FALSE}
#p1
p1_3 
#p2
p2_3 
#p3
p3_3 
#4
p4_3 
#pes1
pes1_3 
#pes2
pes2_3 
```

### Three conditions replication

```{r}
K <- 3
mu <- c(0, 0.4, 0.2)
n <- 82
sd <- 1
design = paste(K, "b", sep = "")
```

```{r, message=FALSE, warning=FALSE}
design_result <- ANOVA_design(
  design = design,
  n = n,
  mu = mu,
  sd = sd,
  labelnames = c("factor1", "level1", "level2", "level3")
  )
```
```{r eval=FALSE, echo=TRUE}
simulation_result <- ANOVA_power(design_result, 
                                 alpha_level = alpha_level, 
                                 nsims = nsims,
                                 verbose = FALSE)
```

```{r echo = FALSE}
knitr::kable(simulation_result_2.15$main_results,
             caption = "Simulated ANOVA Result")%>%
  kable_styling(latex_options = "hold_position")
```


```{r}
exact_result <- ANOVA_exact(design_result,
                            alpha_level = alpha_level,
                            verbose = FALSE)
```


```{r echo=FALSE}
knitr::kable(exact_result$main_results,
             caption = "Exact ANOVA Result")%>%
  kable_styling(latex_options = "hold_position")
```


## Repeated Measures

We can reproduce the same results as Brysbaert finds with his code: 


```{r}
design <- "3w"
n <- 75
mu <- c(0, 0.4, 0.4)
sd <- 1
r <- 0.5
labelnames <- c("speed", "fast", "medium", "slow")
```
We create the within design, and run the simulation

```{r, message=FALSE, warning=FALSE}
design_result <- ANOVA_design(design = design,
                   n = n, 
                   mu = mu, 
                   sd = sd, 
                   r = r, 
                   labelnames = labelnames)
```
```{r eval=FALSE}
simulation_result <- ANOVA_power(design_result, 
                                 alpha_level = alpha_level, 
                                 nsims = nsims,
                                 verbose = FALSE)
```

```{r echo=FALSE}
knitr::kable(simulation_result_3.7$main_results,
             caption = "Simulated ANOVA Result")%>%
  kable_styling(latex_options = "hold_position")
```


```{r}
exact_result <- ANOVA_exact(design_result,
                            alpha_level = alpha_level,
                            verbose = FALSE)
```


```{r echo=FALSE}
knitr::kable(exact_result$main_results,
             caption = "Exact ANOVA Result")%>%
  kable_styling(latex_options = "hold_position")

```

**Results**

The results of the simulation are very similar. Power for the ANOVA *F*-test is around 95.2%. For the three paired t-tests, power is around 92.7. This is in line with the a-priori power analysis when using Gpower:

![](screenshots/gpower_2.png)

We can perform an post-hoc power analysis in Gpower. We can calculate Cohen´s f based on the means and sd, using our own custom formula. 

```{r}
# Our simulation is based onthe following means and sd:
# mu <- c(0, 0.4, 0.4)
# sd <- 1
# Cohen, 1988, formula 8.2.1 and 8.2.2
f <- sqrt(sum((mu - mean(mu)) ^ 2) / length(mu)) / sd 

# We can see why f = 0.5*d.
# Imagine 2 group, mu = 1 and 2
# Grand mean is 1.5, 
# we have sqrt(sum(0.5^2 + 0.5^2)/2), or sqrt(0.5/2), = 0.5.
# For Cohen's d we use the difference, 2-1 = 1. 
```

The Cohen´s *f* is `r f`. We can enter the *f* (using the default 'as in G*Power 3.0' in the option window) and enter a sample size of 75, number of groups as 1, number of measurements as 3, correlation as 0.5. This yields: 

![](screenshots/gpower_3.png)


### Reproducing Brysbaert Variation 1: Changing Correlation

```{r eval=FALSE}
# give sample size
N = 75
# give effect size d
d1 = .4 #difference between the extremes
d2 = .4 #third condition goes with the highest extreme
# give the correlation between the conditions
r = .6 #increased correlation
# give number of simulations
nSim = nsims
# give alpha levels
alpha1 = .05 #alpha level for the omnibus ANOVA
alpha2 = .05 #also adjusted from original by DL
```

```{r eval=FALSE}
# create vectors to store p-values
p1 <- numeric(nSim) #p-value omnibus ANOVA
p2 <- numeric(nSim) #p-value first post hoc test
p3 <- numeric(nSim) #p-value second post hoc test
p4 <- numeric(nSim) #p-value third post hoc test

# define correlation matrix
rho <- cbind(c(1, r, r), c(r, 1, r), c(r, r, 1))
# define participant codes
part <- paste("part",seq(1:N))
for (i in 1:nSim) {
  
  #for each simulated experiment

  data = mvrnorm(n = N,
  mu = c(0, 0, 0),
  Sigma = rho)
  data[, 2] = data[, 2] + d1
  data[, 3] = data[, 3] + d2
  datalong = c(data[, 1], data[, 2], data[, 3])
  conds = factor(rep(letters[24:26], each = N))
  partID = factor(rep(part, times = 3))
  output <- data.frame(partID, conds, datalong)
  test <- aov(datalong ~ conds + Error(partID / conds), 
              data = output)
  tests <- (summary(test))
  p1[i] <- tests$'Error: partID:conds'[[1]]$'Pr(>F)'[[1]]
  p2[i] <- t.test(data[, 1], data[, 2], paired = TRUE)$p.value
  p3[i] <- t.test(data[, 1], data[, 3], paired = TRUE)$p.value
  p4[i] <- t.test(data[, 2], data[, 3], paired = TRUE)$p.value
}
```

```{r eval=FALSE}
sum(p1 < alpha1) / nSim
sum(p2 < alpha2) / nSim
sum(p3 < alpha2) / nSim
sum(p4 < alpha2) / nSim
```
```{r echo = FALSE}
#p1
p1_2
#p2
p2_2
#p3
p3_2
#p4
p4_2
```

```{r}
design <- "3w"
n <- 75
mu <- c(0, 0.4, 0.4)
sd <- 1
r <- 0.6
labelnames <- c("SPEED", 
                "fast", "medium", "slow")
```

We create the 3-level repeated measures design, and run the simulation.

```{r, message=FALSE, warning=FALSE}
design_result <- ANOVA_design(design = design,
                   n = n, 
                   mu = mu, 
                   sd = sd, 
                   r = r, 
                   labelnames = labelnames)
```
```{r eval=FALSE}
simulation_result <- ANOVA_power(design_result, 
                                 alpha_level = alpha_level, 
                                 nsims = nsims,
                                 verbose = FALSE)
```

```{r echo=FALSE}
knitr::kable(simulation_result_3.9$main_results,
             caption = "Simulated ANOVA Result") %>%
  kable_styling(latex_options = "hold_position")
```


```{r}
exact_result <- ANOVA_exact(design_result,
                            alpha_level = alpha_level,
                            verbose = FALSE)
```


```{r echo=FALSE}
knitr::kable(exact_result$main_results,
             caption = "Exact ANOVA Result") %>%
  kable_styling(latex_options = "hold_position")
```



Again, this is similar to GPower for the ANOVA:

![](screenshots/gpower_4.png)