Annotating Bar Plots with Statistical Significance

It only took me a few hours to figure this out. If I'm using 1-way ANOVAs as a statistical test to determine differences in the means of endogenous cortisol concentrations between animals of different sex classes, I can annotate the resulting bar plot with letters denoting significant differences. For example, I'll be using this data set (stored in the variable examp):

   Sample.ID       Sex.Class     EndoF
1          1            Male 134.70700
2          2            Male 112.79506
3          3            Male  83.42866
4          4            Male  66.25100
5          5            Male 111.49403
6          6            Male 105.33600
7          7            Male  68.47000
8          8            Male  68.32500
9          9            Male 100.12880
10        10            Male  90.18200
11        11 Solitary Female  62.52132
12        12 Solitary Female  66.99200
13        13 Solitary Female  77.37900
14        14 Solitary Female  40.61669
15        15 Solitary Female  74.61191
16        16 Solitary Female  29.36257
17        17 Solitary Female  93.74400
18        18 Solitary Female  52.07800
19        19 Solitary Female  50.51210
20        20 Solitary Female  56.12120
21        21    Adult Female  77.55365
22        22    Adult Female  86.24100
23        23    Adult Female 110.45902
24        24    Adult Female  97.06417
25        25    Adult Female 142.45900
26        26    Adult Female 169.85475
27        27    Adult Female 175.39755
28        28    Adult Female 145.52400
29        29    Adult Female  99.12721
30        30    Adult Female 110.12300

I'll calculate and store the means and SEM as per my earlier post. To determine if there are significant differences in the means of EndoF between Sex.Class, I'll punch in the command:

examp.anova = aov(EndoF ~ Sex.Class, data = examp)

This will store the ANOVA data in examp.anova. The summary of the analysis shows:

            Df Sum Sq Mean Sq F value    Pr(>F)  
Sex.Class    2  18666  9333.0  13.502 8.621e-05 ***
Residuals   27  18663   691.2                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

This suggests that there are significant differences! To determine in which sex classes endogenous F are different, I'll run a Tukey pairwise comparison analysis.

examp.pairs <- glht(examp.anova, linfct = mcp(Sex.Class = "Tukey"))

The summary of which shows:

> summary(examp.pairs)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts

Fit: aov(formula = EndoF ~ Sex.Class, data = examp)
Linear Hypotheses:
                                    Estimate Std. Error t value Pr(>|t|)  
Male - Adult Female == 0              -27.27      11.76  -2.319    0.070 .
Solitary Female - Adult Female == 0   -60.99      11.76  -5.187   <0.001 ***
Solitary Female - Male == 0           -33.72      11.76  -2.868    0.021 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

So endogenous F in Solitary Females is significantly different than that in the other classes. I'll plot these results in bar plot form with SEM.

> examp.bp = barplot(examp.means, main = "Endogenous Cortisol by Sex Class", ylab = "Serum Total Cortisol (ng Cortisol/mL Serum)", xlab = "Sex Class", ylim = range(0,200), names.arg = order, axis.lty = 1)

> error.bar(examp.bp,examp.means,examp.sems)

To add letters displaying the statistical significance between groups, I'll just type in the commands:

letters = c("a","a","b") 

text(x=examp.bp,y=examp.means+examp.sems+10,label=letters)

This will add the letters manually to the bar plot 10 units above the SEM whiskers. There's a way to extract the pairwise comparison significance letters automatically using the cld() function.

examp.cld <- cld(examp.pairs)

This will store in examp.cld$mcletters$Letters the appropriate letters. However, they're stored in a different order than in examp.means and examp.sem. You just need to reorder either one. The final product: