3
$\begingroup$

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,

$\endgroup$
0
+50
$\begingroup$

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:

Histogram of density variable, showing large mass at zero and heavy right tail.

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

quantile-quantile plot of the residuals for a linear model on raw density

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

Moving on to answer your question

  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!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.