Solved – Does summary.aov() of MANOVA object adjust for multiple comparisons

anovamanovapost-hoc

I performed a MANOVA analysis in R, testing if two populations differ based on their multivariate phenotype, comprising 16 measurements. I then used the summary.aov() function on the MANOVA model to further assess which measurements differ between the two populations. My question now is, can this ANOVA decomposition serve as a post hoc test and does the summary.aov() function account for multiple testing or should I perform a false discovery rate correction?

To illustrate what I mean, here an example adopted from the R help files:

tear <- c(6.5, 6.2, 5.8, 6.5, 6.5, 6.9, 7.2, 6.9, 6.1, 6.3, 6.7, 6.6, 7.2, 7.1, 6.8, 7.1, 7.0, 7.2, 7.5, 7.6)
gloss <- c(9.5, 9.9, 9.6, 9.6, 9.2, 9.1, 10.0, 9.9, 9.5, 9.4, 9.1, 9.3, 8.3, 8.4, 8.5, 9.2, 8.8, 9.7, 10.1, 9.2)
opacity <- c(4.4, 6.4, 3.0, 4.1, 0.8, 5.7, 2.0, 3.9, 1.9, 5.7, 2.8, 4.1, 3.8,1.6, 3.4, 8.4, 5.2, 6.9, 2.7, 1.9)
Y <- cbind(tear, gloss, opacity)
rate <- factor(gl(2,10), labels = c("Low", "High"))
fit <- manova(Y ~ rate)

This will result in

summary(fit)
          Df  Pillai approx F num Df den Df   Pr(>F)   
rate       1 0.58638    7.561      3     16 0.002273 **
Residuals 18                                           
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

With the ANOVA decomposition being

summary.aov(fit)
Response tear :
            Df Sum Sq Mean Sq F value   Pr(>F)   
rate         1 1.7405 1.74050  12.408 0.002432 **
Residuals   18 2.5250 0.14028                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response gloss :
            Df Sum Sq Mean Sq F value  Pr(>F)  
rate         1 1.3005 1.30050  6.1847 0.02292 *
Residuals   18 3.7850 0.21028                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response opacity :
            Df Sum Sq Mean Sq F value Pr(>F)
rate         1  0.421  0.4205  0.1026 0.7524
Residuals   18 73.785  4.0992

Best Answer

In any situation where you already have a set of unadjusted P values, you can use p.adjust to obtain adjusted ones for any of the methods in p.adjust.methods:

> p.adjust(c(.002432, .02292, .7524), "fdr")
[1] 0.007296 0.034380 0.752400

I had a much more elaborate answer using the lsmeans package, but it does not seem necessary. Let me know if you want to see it.

Related Question