---
title: "Mixed Model vs. within-subjects ANOVA"
date: "`r format(Sys.time(), '%d %B, %Y')`"
author: Henrik Singmann
output: html_notebook
---

The goal of this document is to provide a small simulation comparing repeated-measures ANOVA versus a mixed model in a situation where the maximal-random effects structure for the mixed model is not identified as there is only one observation for each unit of observation (e.g., participant) and cell of the design.

```{r, message=FALSE, results='hide'}
require(afex)
set_sum_contrasts()
require(MASS) 

## functions returning p-value for ANOVA and mixed model
within_anova <- function(data) {
  suppressWarnings(suppressMessages(
  a <- aov_ez(id = "id", dv = "dv", data, within = "condition", return = "univariate", anova_table = list(es = "none"))
  ))
  c(without = a[["univariate.tests"]][2,6],
    gg = a[["pval.adjustments"]][1,2],
    hf = a[["pval.adjustments"]][1,4])
}

within_mixed <- function(data) {
  suppressWarnings(
    m <- mixed(dv~condition+(1|id),data, progress = FALSE)  
  )
  c(mixed=anova(m)$`Pr(>F)`)
}

```

The following code block defines a few parameters for the simulation.

```{r}
NSIM <- 1e4  # number of simulated data sets
NPAR <- 30  # number of participants per cell
NCELLS <- 3  # number of cells (i.e., groups)
```

Next we generate some data without any effect. This is perhaps not done in the most succinct way, but follows the general formulation of a mixed effects model. 

- We first generate an empty `data.frame` containing unit of observation (`id`) and `condition` identifier. 

```{r}
dat <- expand.grid(condition = factor(letters[seq_len(NCELLS)]),
                   id = factor(seq_len(NPAR)))
head(dat)
```
- We create the fixed effects model matrix for our design. As expected, we have one intercept and two fixed effects parameters.

```{r}
mm <- model.matrix(~condition, dat)
head(as.data.frame(mm))
```

- We then define the variance-covariance matrix of the random intercept and random slopes and sample a list of zero-centered individual level parameter from that variance-covariance matrix.
- In the loop we first calculate the dependent variable for each individual participant using the individual random effects estimates estimates (i.e., one row from `d_c_1`).
- Then we add independent residual error to the calculated dv for all data point (of each simulation).

```{r, cache=TRUE}


Sigma_c_1 <- matrix(0.6, NCELLS,NCELLS)
diag(Sigma_c_1) <- 1
d_c_1 <- replicate(NSIM, mvrnorm(NPAR, rep(0, NCELLS), Sigma_c_1), simplify = FALSE)

gen_dat <- vector("list", NSIM)
for(i in seq_len(NSIM)) {
  gen_dat[[i]] <- dat
  gen_dat[[i]]$dv <- NA_real_
  for (j in seq_len(NPAR)) {
    gen_dat[[i]][(j-1)*3+(1:3),"dv"] <- mm[1:3,] %*% d_c_1[[i]][j,]
  }
  gen_dat[[i]]$dv <- gen_dat[[i]]$dv+rnorm(nrow(mm), 0, 1)
}

```

After generating the data, we only need to apply the ANOVA and mixed model to each generated data set.

```{r, cache=TRUE}
p_c1_within <- vapply(gen_dat, within_anova, rep(0.0, 3))
m_c1_within <- vapply(gen_dat, within_mixed, 0.0)
```

The histogram of the p-value distribution clearly shows that the mixed model and the standard ANOVA have basically the same p-value distribution.

```{r, fig.width=7, fig.height=5}
ylim <- c(0, 700)
par(mfrow = c(1,3))
hist(p_c1_within[1,], breaks = 20, main = "ANOVA (default)", xlab = "p-value", ylim=ylim)
hist(p_c1_within[2,], breaks = 20, main = "ANOVA (GG)", xlab = "p-value", ylim=ylim)
hist(m_c1_within, breaks = 20, main = "Random-Intercept Model", xlab = "p-value", ylim=ylim)

```

This is also visible when looking at the proportion of p-values smaller than .05:

```{r}
mean(p_c1_within[1,] < 0.05)
mean(p_c1_within[2,] < 0.05)
mean(p_c1_within[3,] < 0.05)
mean(m_c1_within < 0.05)
```
