# Running a simple anova test in R (repeated measures)

We are studying 3 different proteins, each under 9 different conditions at 3 different timepoints (Day 1,2,3). For these we have 3 biological replicates. So we have 81 different experiments (3 proteins * 9 conditions * 3 replicates) ‒ and for each experiment we have data at three different timepoint readings on consecutive days. This gives us 243 observations in a balanced design.

We would like to show which of these proteins and conditions are statistically different from each other. We would like a comparision between proteins, and the conditions of each protein compared. For this we were thinking of using a repeated measure anova test (using R).

I replicated a MWE of the dataset and example here:

library(RCurl)
library(dplyr)

raw.data <- getURL("https://gist.githubusercontent.com/jp-um/1849ac4ac61411d0751cdbec4406e0cd/raw/4b014f986085665e75806c38a25f39093b2d19df/anon.csv")
exp.data <- read.csv(text = raw.data, colClasses=c("experiment"="factor",
"protein"="factor",
"condition"="factor",
"day"="factor",
"bioreplicate"="factor"))
summary(exp.data, maxsum=10)

aov.model <- aov(density ~ protein*condition*day + Error(experiment/day), data = exp.data)
summary(aov.model)


The output is:

Error: experiment
Df   Sum Sq Mean Sq F value Pr(>F)
protein            2 10729989 5364994  1166.3 <2e-16 ***
condition          8 16430568 2053821   446.5 <2e-16 ***
protein:condition 16 29649758 1853110   402.8 <2e-16 ***
Residuals         54   248404    4600
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: experiment:day
Df  Sum Sq Mean Sq F value  Pr(>F)
day                     2   92976   46488   13.24 7.2e-06 ***
protein:day             4 1776592  444148  126.49 < 2e-16 ***
condition:day          16 3419459  213716   60.87 < 2e-16 ***
protein:condition:day  32 7415908  231747   66.00 < 2e-16 ***
Residuals             108  379221    3511
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


I have a few questions, please:

1. Is repeated measure anova the way to go here (as opposed to mixed models)?
2. Is this the correct way to specify the formula? Will this give me repeated measures over the timepoints (Day)? Specifically what does Error(experiment/day) mean? My interpretation is that we have a random effect based on the biological replicate (experiment in my case) and repeated readings for the timepoint (Day).
3. What is the ezANOVA equivalent way to write this?
4. Is the above anova, equivalent to the linear mixed effect model lme(density ~ protein*condition*day, random = ~1 | experiment/day, data = exp.data) ?
5. The output tells me that there are differences between them, but I would like to know which combination gives the difference. I know I can use a post-hoc test for this, and I found TukeyHSD does not work on repeated measures. I have found I can use glht for this; but I am unable to interpret its output. (I tried glht(lme.model, linfct=mcp(protein="Tukey", condition="Tukey")) but I am not sure this is correct)

Apologies for the (many, and rather basic I am afraid) questions. I would appreciate your time and help.

Many thanks,

## Before I begin, the most important thing for you to do is log-transform density in your model.

That's because density is bounded at zero, and has a very heavy right tail:

Furthermore, the Q-Q plot of your model residuals reveals very heavy tails:

This is an absolutely clear example of needing a log transformation. It happens to resolve these problems in your case.

1. Repeated measures ANOVA is fine here (so would be mixed-effect modeling).
2. This is the way. You would get the same result if you rewrote the error term as Error( experiment ). This tells aov() to partition the residual error into two parts for the analysis of variance: differences between experiments, and differences within experiments. The "within experiments" part is the sum of squared differences between the observed and fitted values for each day's density. The "between experiments" part is the sum of squared differences between the observed vs fitted mean density of the experiments.
3. The equivalent ezANOVA() call is:
ezANOVA( exp.data3, dv= "density", wid=experiment, within=day, between=c("protein", "condition") )

1. Your aov() and lme() models are not exactly the same, but they are very similar. They're so similar because the design is balanced and the variance is quite equal between experiments. If the design was unbalanced or the variance was not so equal, then the partial pooling of residual error by lme() would become more important.
2. In order to get results from glht() you need to assign its result to a variable, and then summarize that variable:
tmp = glht(lme.model, linfct=mcp(protein="Tukey", condition="Tukey"))
summary( tmp )

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: lme.formula(fixed = log(density) ~ protein * condition * day,
data = exp.data, random = ~1 | experiment)

Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
protein: DEF - ABC == 0  0.80355    0.12817   6.269    <0.01 ***
protein: GHI - ABC == 0 -1.32768    0.12817 -10.358    <0.01 ***
protein: GHI - DEF == 0 -2.13122    0.12817 -16.628    <0.01 ***
condition: lm - jk == 0  0.13706    0.12817   1.069   0.9878
condition: no - jk == 0  0.11786    0.12817   0.920   0.9962
condition: pq - jk == 0  0.16504    0.12817   1.288   0.9555
condition: rs - jk == 0  1.29121    0.12817  10.074    <0.01 ***
condition: tu - jk == 0  0.07276    0.12817   0.568   0.9999
condition: vw - jk == 0 -0.23187    0.12817  -1.809   0.7167
condition: xy - jk == 0 -0.28587    0.12817  -2.230   0.4164
condition: z_ - jk == 0 -0.38232    0.12817  -2.983   0.0772 .
condition: no - lm == 0 -0.01920    0.12817  -0.150   1.0000
condition: pq - lm == 0  0.02798    0.12817   0.218   1.0000
(etc)


Now it is easier to see the comparisons being made, as well as the standard error for the comparisons and the resulting z-statistic and p-value. That said, I prefer to use emmeans() over glht() because to me it is more explicit about what values are used for the other factors in the comparison. But I am not totally sure how their computations differ (I haven't dug that deeply into the code, myself).

> emmeans(aov.model,  pairwise ~ protein + condition)

Note: re-fitting model with sum-to-zero contrasts
NOTE: Results may be misleading due to involvement in interactions
<snip>
\$contrasts
contrast        estimate     SE df t.ratio p.value
ABC jk - DEF jk -1.29127 0.0945 54 -13.668  <.0001
ABC jk - GHI jk  0.80242 0.0945 54   8.494  <.0001
ABC jk - ABC lm -0.37632 0.0945 54  -3.983  0.0395
ABC jk - DEF lm -3.02000 0.0945 54 -31.967  <.0001
ABC jk - GHI lm -1.98646 0.0945 54 -21.027  <.0001
ABC jk - ABC no  0.06069 0.0945 54   0.642  1.0000
ABC jk - DEF no -1.12982 0.0945 54 -11.959  <.0001
ABC jk - GHI no  0.05284 0.0945 54   0.559  1.0000
ABC jk - ABC pq -0.09484 0.0945 54  -1.004  1.0000
ABC jk - DEF pq -1.52451 0.0945 54 -16.137  <.0001
ABC jk - GHI pq -2.47809 0.0945 54 -26.231  <.0001
<snip>
Results are averaged over the levels of: day
Results are given on the log (not the response) scale.
P value adjustment: tukey method for comparing a family of 27 estimates


I hope this helps!