[latex]\newcommand{\pr}[1]{P(#1)} \newcommand{\var}[1]{\mbox{var}(#1)} \newcommand{\mean}[1]{\mbox{E}(#1)} \newcommand{\sd}[1]{\mbox{sd}(#1)} \newcommand{\Binomial}[3]{#1 \sim \mbox{Binomial}(#2,#3)} \newcommand{\Student}[2]{#1 \sim \mbox{Student}(#2)} \newcommand{\Normal}[3]{#1 \sim \mbox{Normal}(#2,#3)} \newcommand{\Poisson}[2]{#1 \sim \mbox{Poisson}(#2)} \newcommand{\se}[1]{\mbox{se}(#1)} \newcommand{\prbig}[1]{P\left(#1\right)} \newcommand{\degc}{$^{\circ}$C}[/latex]

The one-way analysis of variance [latex]F[/latex] test can be used to show that there is a difference between the mean responses from multiple groups. For example, in Chapter 19 we saw very strong evidence that the mean water loss was different for the three wind speeds in the water loss data table. The obvious question is to then ask where the differences lie.

Least Significant Difference

We will start with an example which has a moderate number of possible differences. A dataset in the Appendix gives the yields from six treatments applied to potatoes in an agricultural trial. The treatments themselves are combinations of factors and were arranged in a latin square but for now we will simply look for differences in the mean yield between the six treatments.

The figure above shows the distribution of potato yield for each treatment group. Here the null hypothesis is
\[ H_0: \mu_A = \mu_B = \mu_C = \mu_D = \mu_E = \mu_F, \]
that the six groups have equal mean yield. We can carry out a one-way analysis of variance to test this, giving the ANOVA table in below.

ANOVA table for potato yield (lbs) by treatment

From the [latex]F[/latex] distribution table, the critical value for significance at the 0.001 level is [latex]F^{*}_{5,30} = 5.53[/latex] so this is strong evidence to suggest that the means are not all equal. But which ones are different?

One approach to answering this is to use a two-sample [latex]t[/latex] test to compare differences between each pair of treatments. For example, we could test the null hypothesis [latex]H_0: \mu_A = \mu_B[/latex] by comparing the means [latex]\overline{x}_A[/latex] and [latex]\overline{x}_B[/latex], as given in the following table. Since analysis of variance assumes common standard deviations, we can use a pooled two-sample [latex]t[/latex] test and, moreover, we can use the pooled estimate of variability from all six treatments since this gives more information about this common standard deviation than just combining the variability estimates from [latex]A[/latex] and [latex]B[/latex]. The pooled variance is [latex]\mbox{MSR} = 3640[/latex] so the pooled estimate of standard deviation is
\[ \sqrt{\mbox{MSR}} = \sqrt{3640} = 60.33 \mbox{ lbs}, \]
which seems reasonable given the individual standard deviations in the table below.

Summary statistics for potato yield (lbs) by treatment

Our [latex]t[/latex] statistic for comparing [latex]A[/latex] and [latex]B[/latex] is thus
\[ t_{30} = \frac{(345.0 – 426.5) – 0}{60.33\sqrt{\frac{1}{6} + \frac{1}{6}}} = -2.340. \]
Note again that the degrees of freedom of this statistic are 30 instead of 10 since we are pooling together variability across the six treatments. From Student’s T distribution we find that the one-sided [latex]P[/latex]-value is between 0.025 and 0.01, giving a two-sided [latex]P[/latex]-value between 0.05 and 0.02. Thus there is some evidence of a difference in mean potato yield between treatments [latex]A[/latex] and [latex]B[/latex].

We could continue doing [latex]t[/latex] tests for each of the other pairs of treatments. In general, for [latex]m[/latex] groups there are
\[{m \choose 2} = \frac{m(m-1)}{2} \]
tests to try. With [latex]m = 6[/latex] treatments this gives 15 pairs, as shown in the table below. This would be straightforward but certainly tedious.

Pairwise differences in mean yield (lbs)

However we can simplify this task by noting that the denominator of the [latex]t[/latex] statistic above will be the same for each test since we are using the common pooled estimate and since we have the same number of observations in each group. The only way we can get two-sided significance at the 5% level is if the difference in means divided by this number is at least the critical [latex]t[/latex] value, [latex]t_{30}^{*} = 2.042[/latex]. That is, we require
\[ \frac{\mbox{difference}}{60.33\sqrt{\frac{1}{6} + \frac{1}{6}}} \ge 2.042. \]
What is the smallest difference that would be significant? Rearranging this equation gives
\[ \mbox{difference} \ge 2.042 \times 60.33\sqrt{\frac{1}{6} + \frac{1}{6}} = 71.1 \mbox{ lbs}. \]
This value is called the Least Significant Difference (LSD) and gives a simple tool for assessing all the pairwise comparisons. We’ve already seen that the difference between [latex]A[/latex] and [latex]B[/latex] was 81.5 lbs; this is more than the LSD and so there is evidence at the 5% level of a difference between the corresponding means, as noted before. Now we can look through all the remaining differences in the previous table and decide which are significant. We find all differences to be significant at the 5% level except for [latex]A-D[/latex], [latex]B-C[/latex], [latex]B-D[/latex] and [latex]C-E[/latex]. The difference [latex]C-D[/latex] is only just significant.

We could likewise calculate a 95% confidence interval for each difference since the LSD is precisely the margin of error we would use. For example, a 95% confidence interval for [latex]A-B[/latex] is
\[ -81.5 \pm 71.1, \]
or [latex](-152.6,-10.4)[/latex]. Looking to see if 0 is in this range is equivalent to comparing the mean difference to the LSD.

Although the LSD simplifies the testing process, there is still a problem with making multiple comparisons in this way. Recall that testing for significance at the 5% level means that there is a 5% chance of making a mistake, rejecting the null hypothesis when there is no real effect (a Type I error). Correspondingly there is a 95% chance of not making a mistake. With 15 tests the probability of making no mistakes is [latex]0.95 \times \cdots \times 0.95 = 0.95^{15} = 0.4633[/latex] (supposing the tests are independent). So the chance of making at least one mistake is now 53.67%, much higher than the specified level of 5%.

Note that there is some safety in doing these pairwise comparisons only after having found evidence of an overall difference from the analysis of variance. That is, if we had not found evidence from the [latex]F[/latex] test then we would not have proceeded to do all these comparisons and so we would have been protected from the possibility of making these Type I errors. This two-stage approach is known as the protected Fisher’s least significant difference test. Hayter (1986) showed that by doing this the overall error rate was no more than around 29% when comparing six groups, better than 53.67% but still very high.

Beyond the protected Fisher’s LSD, there are several alternative ways of correcting for this issue in multiple comparisons (Glover & Mitchell, 2002; Kleinbaum, Kupper, & Muller, 1988; Steel & Torrie, 1960). We will describe two which give some of the basic ideas involved.

Bonferroni Correction

A simple approach is to adjust the pairwise error rate so that the family error rate is no more than some target level. For example, if we want the family error rate to be no more than 0.05 and we are carrying out [latex]k[/latex] pairwise tests, a Bonferroni correction says to look for significance at the
\[ \alpha = \frac{0.05}{k} \]
level in each test. In this case, 0.05/15 = 0.0033, so we would only reject an individual null hypothesis if the [latex]P[/latex]-value was less than 0.0033. Similarly, to find confidence intervals with a family error rate of 5% we could calculate each individual interval at 99.67% confidence instead of 95%.

The following table shows the [latex]P[/latex]-values obtained from software for the two-sample [latex]t[/latex] tests for each pair of treatments, as in the previous section. Although the [latex]P[/latex]-values for [latex]A - B[/latex], [latex]B - E[/latex], [latex]C - D[/latex] and [latex]E - F[/latex] are all below 5%, using the Bonferroni correction we would not take any of these comparisons as significant. Other comparisons, such as [latex]D - E[/latex], give [latex]P[/latex]-values below 0.0033 and thus are kept as significant differences.

Pairwise [asciimath]t[/asciimath] tests for differences in mean yield (lbs)

The critical value for 99.67% confidence is [latex]t_{30}^{*} = 3.189[/latex], compared to 2.042 for 95% confidence, and we can use this to calculate Bonferroni confidence intervals for each difference.
For example, a 99.67% confidence interval for [latex]A-B[/latex] is
\[ -81.5 \pm 3.189 \times 60.33\sqrt{\frac{1}{6} + \frac{1}{6}} = -81.5 \pm 111.1, \]
or [latex](-192.6, 29.6)[/latex]. Looking to see if 0 is in this range is equivalent to finding [latex]P[/latex]-values below 0.0033. Clearly 111.1 lbs is much larger than the LSD of 71.1 lbs in the previous section, reflecting the adjustment for multiple comparisons. The [latex]A-B[/latex] difference was significant for the LSD test but is no longer significant after the Bonferroni correction.

The Bonferroni method is very simple but is quite conservative, particularly for larger numbers of comparisons.
Note that this is also the worst-case scenario in the sense that it is based on making all possible comparisons. Frequently we are interested in any significant difference, but sometimes we may have a priori reasons to only be looking at selected differences. For example, in this experiment we might have taken treatment [latex]A[/latex] as a baseline and wanted to compare the other five treatments with it. Thus there would only be 5 comparisons to consider, instead of 15, and so the Bonferroni correction would only require significance at the [latex]0.05/5 = 0.01[/latex] level.

Microarray Data

Multiple comparison procedures are very important in the analysis of genetic microarray data (Leung & Cavalieri, 2003). A common aim of microarrays is to identify differentially expressed genes between groups and they allow thousands of genes to be compared simultaneously. The immediate problem then is that we would expect a large number of spurious differences to be detected.

The Islanders have a shorter genome than regular humans and so make a manageable example to illustrate the principles of microarray analysis. Chloe Moore from Colmar University was interested in identifying any Islander genes related to the risk of diabetes. She recruited 20 subjects for her study, 10 without diabetes and 10 with diabetes. For each subject she conducted a microarray experiment. The figure below shows the resulting expression levels for each of the 10 subjects with diabetes and for each of the 10 subjects without diabetes. The question is whether there is any systematic difference between the intensity levels expressed in each set of images.

The figure below shows a more quantitative summary of the results with side-by-side box plots for each of the 256 genes appearing in the microarrays. For clarity we have indicated the median location with a horizontal line, rather than the dot we normally use. Look around the box plots and see if you can identify any genes where there may be a significant difference between the expression levels. (We will label the genes starting with B000 at the top left, going across to B015 at the top right and then proceeding down the plot to B255 in the bottom right.)

For each of the 256 genes we can carry out a two-sample [latex]t[/latex] test to compare the levels between the subjects with diabetes and those without. Each test gives a [latex]P[/latex]-value and these can be visualised with a Manhattan plot, seen below. In this plot the horizontal axis gives the gene position while the vertical axis gives the value [latex]-\log_{10}(p)[/latex] for each [latex]P[/latex]-value, [latex]p[/latex]. This transformation means that smaller [latex]P[/latex]-values will appear higher on the scale. For example, the highest value of around 3.0 corresponds to a [latex]P[/latex]-value of around [latex]10^{-3.0} = 0.001[/latex].

Four genes with strongest [asciimath]P[/asciimath]-values

Gene	[asciimath]p[/asciimath]
B116	0.000998
B041	0.001415
B186	0.001559
B118	0.005669

The table above shows the genes with the four strongest [latex]P[/latex]-values in terms of the difference in mean expression levels between the two groups and the following figure shows the corresponding box plots in more detail. The strongest of these, B116, has a very small [latex]P[/latex]-value of 0.000998. As an individual result we would take this to be substantial evidence that the expression of this gene was related to diabetes, with diabetic subjects having lower levels of expression than non-diabetic subjects. However, we have carried out 256 hypothesis tests — would we expect to find a [latex]P[/latex]-value this small by chance?

In our discussion of the Bonferroni correction above we were taking a fixed error rate and adjusting it to find the new error rate we should use in our analysis. Rather than compare [latex]P[/latex]-values to an adjusted [latex]\alpha[/latex] value, an alternative is to adjust the [latex]P[/latex]-values themselves (Wright, 1992). Here instead of dividing a target [latex]\alpha[/latex] by the number of tests, [latex]k = 256[/latex], we multiply the individual [latex]P[/latex]-values by 256 to obtain the adjusted P-value. The table below gives the results of this simple calculation and we see that the evidence is now much weaker. In particular, evidence for a B116 effect is now inconclusive.

Four genes with strongest adjusted [asciimath]P[/asciimath]-values

Gene	[asciimath]p[/asciimath] (adjusted)
B116	0.255
B041	0.362
B186	0.399
B118	1.000

Of course this adjustment is very conservative. All but three of the genes had [latex]P[/latex]-values that when multiplied by 256 gave probabilities above 1 (and so have been capped at 1). Wright (1992) discusses more sophisticated adjustments to multiple [latex]P[/latex]-values that you will find available in statistical software though the basic principle of correcting for a family-wise error rate is still the same. Reiner et al. (2003) gives an overview of alternative methods based on adjusting for the false discovery rate, the expected probability of false positives when looking for factors affecting gene expression.

Note that even though the results of this microarray study have been inconclusive they do help direct further research. B116 and the other strongest genes can be taken as candidates for more detailed investigation.

Tukey’s Honestly Significant Difference

The problem with making conclusions from a posteriori comparisons, those that we think of only after looking at the data, is that we will look for the largest difference between sample means. With lots of pairwise comparisons it can be quite likely this largest difference is significant even if all the population means were the same. The Bonferroni correction is a simple way of dealing with this, reducing the individual error rate to compensate. A more sophisticated approach is given by Tukey (1953).

The idea behind Tukey’s method is to actually use the sampling distribution of the largest standardised difference found when comparing [latex]k[/latex] groups which all have the same means and standard deviations. This distribution is called the Studentised range distribution, [latex]Q(k,d)[/latex], where [latex]k[/latex] is the number of groups and [latex]d[/latex] is the residual degrees of freedom. The method described here is only valid for the case of equal replications in each group. Steel and Torrie (1960) give modifications for the more general case.

The Tukey statistic is a minor modification of the regular [latex]t[/latex] statistic with
\[ q_{k,d} = \frac{|\overline{x}_i – \overline{x}_j|}{\sqrt{\frac{\mbox{MSR}}{n}}}, \]
where [latex]n[/latex] is the number of replications in each group. For example, for comparing treatments [latex]C[/latex] and [latex]F[/latex] in the potato experiment we would have
\[ q_{6,30} = \frac{124}{60.33/\sqrt{6}} = 5.03. \]
Instead of comparing to the [latex]t[/latex] distribution, we obtain bounds on the [latex]P[/latex]-value from the critical values of the Studentised range distribution at the end of this chapter. Here we find that the [latex]P[/latex]-value is between 0.05 and 0.01, some evidence of a difference between [latex]C[/latex] and [latex]F[/latex] in their mean potato yield.

As with the two-sample [latex]t[/latex] test, the denominator of this statistic is the same for each pairwise comparison and so we can determine the minimum difference that would be significant. The critical value for 5% significance from the Studentised range distribution is [latex]q_{6,30}^{*} = 4.30[/latex] so that the difference must be at least
\[ 4.30 \times \frac{60.33}{\sqrt{6}} = 105.9 \mbox{ lbs}. \]
This value is known as Tukey’s Honestly Significant Difference (HSD). Note that is it higher than the least significant difference, 71.1 lbs, found earlier in this chapter, taking into account the multiple comparisons by requiring a higher threshold for significance. However it is less conservative than the Bonferroni correction where the corresponding value was 111.1 lbs.

You can refer back to the list of pairwise differences in the earlier table and compare each one to the HSD to assess significance. Alternatively, the table below used software to calculate the [latex]P[/latex]-values from the Studentised range distribution for each pairwise difference.

Pairwise Tukey tests for differences in mean yield (lbs)

Note that if we were just comparing two groups then the critical value would be [latex]q_{2,30}^{*} = 2.89[/latex] and so the HSD would be
\[ 2.89 \times \frac{60.33}{\sqrt{6}} = 71.1 \mbox{ lbs}. \]
This is just the LSD, again illustrating that Tukey’s HSD is taking into account the multiple comparisons being performed by increasing the threshold for significance.

The advantage of the LSD and Bonferroni methods is that they can be performed using standard [latex]t[/latex] procedures. The Bonferroni method may also be preferable if there is a valid a priori limit on the number of comparisons being made. However, in software the Tukey procedure is often the standard approach available.

Studentized Range distribution [latex]Q(k,d)[/latex]

This table gives gives [latex]q^{*}[/latex] such that [latex]\pr{Q_{k,d} \ge q^{*}} = p[/latex].