Multiple Comparisons in Analysis of Variance
StatsDirect provides functions for multiple comparison (simultaneous inference), specifically all pairwise comparisons and all comparisons with a control. For k groups there are k(k-1)/2 possible pairwise comparisons.
Tukey (Tukey-Kramer if unequal group sizes), Scheffé, Bonferroni and Newman-Keuls methods are provided for all pairwise comparisons (Armitage and Berry, 1994; Wallenstein, 1980; Miller, 1981; Hsu, 1996; Kleinbaum et al., 1998). Dunnett's method is used for multiple comparisons with a control group (Hsu, 1996).
For k groups, ANOVA can be used to look for a difference across k group means as a whole. If there is a statistically significant difference across k means then a multiple comparison method can be used to look for specific differences between pairs of groups. The reason that two sample methods should not be used to make multiple pairwise comparisons is that they are not designed for repeat testing in a "data dredging" manner.
If 20 repeat pairwise tests are made then you can not accept the conventional 1 in 20 chance of being wrong as a cut off level for statistical inference, i.e. there is a higher risk of type I error. A simple solution to this problem is to reduce the cut-off for statistical significance with increasing numbers of contrasts made; Bonferroni's method does just this with multiple t tests. More sophisticated methods, such as Tukey(-Kramer), consider the statistical distributions associated with systematic repeated testing; both Tukey(-Kramer) and Newman-Keuls methods are based upon the Studentized range statistic. Scheffé's method gives a very conservative/cautious weighting against the risk of type I error and is therefore less powerful for the detection of true differences. The most acceptable general method for all pairwise comparisons is Tukey(-Kramer), the P values for which are exact with balanced designs (Hsu, 1996).
The contrasts in each table are listed in decreasing order of the absolute standardized difference between the two group means, |L/SE(L)|; with equal group sizes this is the order of the differences themselves, and the Newman-Keuls and Dunnett tables are ordered by the difference. In the Tukey and Scheffé tables the word "stop" is shown next to the first non-significant P value: you should not consider further contrasts if you are making a simultaneous analysis (similar to the Shaffer-Holm method). In the Newman-Keuls table each P value is referred to the range of ordered means that its pair spans, so the P values are not in order; a pair that lies inside a wider range that is not significant is marked "not tested" and does not count as different whatever its own P value, and the summary of significant contrasts shows which pairs differ.
The following is a decision tree for selecting a multiple contrast method:
- pairwise
- equal groups sizes: Tukey
- unequal group sizes: Tukey-Kramer or Scheffé
- not pairwise
- with a control: Dunnett
- planned: Bonferroni
- not planned: Scheffé
Note that Bonferroni and Scheffé methods are completely general; they can be used for unplanned (a posteriori) or planned (a priori) multiple comparisons.
This is a controversial area in statistics and you would be wise to seek the advice of a statistician at the design stage of your study. In general you should design experiments so that you can avoid having to "dredge" groups of data for differences, decide which contrasts you are interested in at the outset. Note that multiple independent comparisons (e.g. multiple t or Mann-Whitney tests) may be justified if you identify the comparisons as valid at the design stage of your investigation.
Other statistical software may refer to LSD (least significant difference) methods, please note that the Bonferroni technique described above is an LSD method.
Example
Test workbook (ANOVA worksheet: Substance 1, Substance 2, Substance 3, Substance 4).
| Substance 1 | Substance 2 | Substance 3 | Substance 4 |
| 29 | 17 | 17 | 18 |
| 28 | 25 | 16 | 20 |
| 23 | 24 | 21 | 25 |
| 26 | 19 | 22 | 24 |
| 26 | 28 | 23 | 16 |
| 19 | 21 | 18 | 20 |
| 25 | 20 | 20 | 20 |
| 29 | 25 | 17 | 17 |
| 26 | 19 | 25 | 19 |
| 28 | 24 | 21 | 17 |
This illustration uses four columns of the test workbook, made up for the purpose: a response measured ten times under each of four substances, as tabulated above. Select One Way from the Analysis of Variance section of the analysis menu, select the four columns in one action, and then choose each multiple comparison method in turn from the options offered after the analysis (Bonferroni asks for the two groups to compare, here Substance 1 and Substance 4, and for the number of comparisons planned, here 6 for all the pairs; Dunnett asks for the control group, here Substance 1).
For this example:
One way analysis of variance
Variables: Substance 1, Substance 2, Substance 3, Substance 4
| Source of Variation | Sum Squares | DF | Mean Square |
| Between Groups | 249.875 | 3 | 83.291667 |
| Within Groups | 350.9 | 36 | 9.747222 |
| Corrected Total | 600.775 | 39 |
F (variance ratio) = 8.54517 P = 0.0002
Tukey multiple comparisons
Critical value (Studentized range) = 3.808798, |q*| = 2.693227
Pooled standard deviation = 3.122054
| Comparison | Mean difference L (95% CI) | q | |
| Substance 1 vs. Substance 4 | 6.3 (2.539649 to 10.060351) | 6.381167 | P = 0.0004 |
| Substance 1 vs. Substance 3 | 5.9 (2.139649 to 9.660351) | 5.976014 | P = 0.0009 |
| Substance 1 vs. Substance 2 | 3.7 (-0.060351 to 7.460351) | 3.74767 | P = 0.0552 stop |
| Substance 2 vs. Substance 4 | 2.6 (-1.160351 to 6.360351) | 2.633498 | P = 0.2621 |
| Substance 2 vs. Substance 3 | 2.2 (-1.560351 to 5.960351) | 2.228344 | P = 0.4049 |
| Substance 3 vs. Substance 4 | 0.4 (-3.360351 to 4.160351) | 0.405153 | P = 0.9917 |
| Variable | Mean | Significant contrasts |
| Substance 1 | 25.9 | Substance 4, Substance 3 |
| Substance 4 | 19.6 | Substance 1 |
| Substance 3 | 20 | Substance 1 |
| Substance 2 | 22.2 | none |
The statistic tabulated for Tukey and Newman-Keuls is the Studentized range q: the difference between the two means divided by the standard error of a mean, s/√n (or its Tukey-Kramer form for unequal sizes), which is compared with the critical value. Scheffé's statistic is the difference divided by the standard error of the difference, s√(2/n), compared with its own critical value.
Scheffé multiple comparisons
Critical value = 2.93237
| Comparison | Mean difference L (95% CI) | |L/SE(L)| | |
| Substance 1 vs. Substance 4 | 6.3 (2.205751 to 10.394249) | 4.512167 | P = 0.001 |
| Substance 1 vs. Substance 3 | 5.9 (1.805751 to 9.994249) | 4.22568 | P = 0.0021 |
| Substance 1 vs. Substance 2 | 3.7 (-0.394249 to 7.794249) | 2.650003 | P = 0.0896 stop |
| Substance 2 vs. Substance 4 | 2.6 (-1.494249 to 6.694249) | 1.862164 | P = 0.34 |
| Substance 2 vs. Substance 3 | 2.2 (-1.894249 to 6.294249) | 1.575677 | P = 0.4874 |
| Substance 3 vs. Substance 4 | 0.4 (-3.694249 to 4.494249) | 0.286487 | P = 0.9938 |
| Variable | Mean | Significant contrasts |
| Substance 1 | 25.9 | Substance 4, Substance 3 |
| Substance 4 | 19.6 | Substance 1 |
| Substance 3 | 20 | Substance 1 |
| Substance 2 | 22.2 | none |
Newman-Keuls multiple comparisons
| Comparison | Mean difference L | Separation | q | |
| Substance 1 vs. Substance 4 | 6.3 | 4 | 6.381167 | P = 0.0004 |
| Substance 1 vs. Substance 3 | 5.9 | 3 | 5.976014 | P = 0.0004 |
| Substance 1 vs. Substance 2 | 3.7 | 2 | 3.74767 | P = 0.0119 |
| Substance 2 vs. Substance 4 | 2.6 | 3 | 2.633498 | P = 0.1644 |
| Substance 2 vs. Substance 3 | 2.2 | 2 | 2.228344 | P = 0.1238 not tested |
| Substance 3 vs. Substance 4 | 0.4 | 2 | 0.405153 | P = 0.7761 not tested |
| Variable | Mean | Significant contrasts |
| Substance 1 | 25.9 | Substance 4, Substance 3, Substance 2 |
| Substance 4 | 19.6 | Substance 1 |
| Substance 3 | 20 | Substance 1 |
| Substance 2 | 22.2 | Substance 1 |
Bonferroni comparison
Variable A: Substance 1
Variable B: Substance 4
Mean A - Mean B = 6.3
Estimated std. error = 1.396225
Number of groups = 4
95% confidence interval = 3.468324 to 9.131676
Bonferroni-adjusted (simultaneous) 95% confidence interval for 6 comparisons = 2.401779 to 10.198221 (each comparison at 99.17%)
t = 4.512167
df = 36
P < 0.0001
Bonferroni critical P for 6 comparisons = 0.008333
Dunnett multiple comparisons with a control
Critical value (|d|) = 2.452127
Pooled standard deviation = 3.122054
Control (n) = Substance 1 (10)
| Comparison (n) | Mean difference (95% CI) | |
| Substance 4 (10) | -6.3 (-9.723721 to -2.876279) | P = 0.0002 |
| Substance 3 (10) | -5.9 (-9.323721 to -2.476279) | P = 0.0004 |
| Substance 2 (10) | -3.7 (-7.123721 to -0.276279) | P = 0.0316 |
Every method finds the response to Substance 1 higher than to Substances 3 and 4. Whether it also differs from Substance 2 depends on the method: Newman-Keuls says so (P = 0.0119), Tukey (P = 0.0552) and Scheffé (P = 0.0896) do not, and the word "stop" marks where each of those two sequences of contrasts ends; in the Newman-Keuls table the two pairs inside the non-significant range from Substance 4 to Substance 2 are marked as not tested. With Substance 1 as the control, Dunnett's method finds all three of the other substances to give a lower response.
R code
This R code reproduces the example above. It needs no packages and was checked with R 4.6.1. Paste it into R, or save it as a script and run it.
# Multiple comparisons after one way analysis of variance: an illustration with
# the Substance 1 to 4 columns of the StatsDirect test workbook (ANOVA worksheet)
s1 <- c(29, 28, 23, 26, 26, 19, 25, 29, 26, 28)
s2 <- c(17, 25, 24, 19, 28, 21, 20, 25, 19, 24)
s3 <- c(17, 16, 21, 22, 23, 18, 20, 17, 25, 21)
s4 <- c(18, 20, 25, 24, 16, 20, 20, 17, 19, 17)
y <- c(s1, s2, s3, s4)
group <- factor(rep(paste("Substance", 1:4), each = 10))
fit <- aov(y ~ group)
print(summary(fit))
# Tukey: R's TukeyHSD gives every pairwise difference with its simultaneous
# confidence interval and P value, the figures in StatsDirect's table, but with
# the opposite sign (R subtracts the first level from the second, so the limits
# are swapped too) and in a different order
print(TukeyHSD(fit))
# The pieces of the table, calculated. The groups are of equal size here; the
# Tukey-Kramer form for unequal sizes replaces 1 / n by (1 / ni + 1 / nj) / 2
# in the standard error
six <- function(x) formatC(x, digits = 6, format = "f", drop0trailing = TRUE)
pv <- function(p) {
if (p < 0.0001) "P < 0.0001" else
paste("P =", formatC(p, digits = 4, format = "f", drop0trailing = TRUE))
}
a <- anova(fit)
mse <- a$"Mean Sq"[2] # the pooled variance
nu <- a$Df[2] # its degrees of freedom, 36
k <- nlevels(group)
n <- 10
means <- tapply(y, group, mean)
pairs <- combn(k, 2) # every pair i < j, 6 of them
delta <- means[pairs[1, ]] - means[pairs[2, ]]
ord <- order(-abs(delta)) # the report's order
label <- paste(names(means)[pairs[1, ]], "vs.", names(means)[pairs[2, ]])
q <- qtukey(0.95, k, nu)
cat("Tukey: Critical value (Studentized range) =", six(q), " |q*| =", six(q / sqrt(2)),
"\n")
cat("Pooled standard deviation =", six(sqrt(mse)), "\n")
se <- sqrt(mse / n) # the difference divided by this is the
for (i in ord) { # Studentized range statistic
cat(label[i], six(delta[i]), six(delta[i] - q * se), "to", six(delta[i] + q * se),
six(abs(delta[i]) / se),
pv(ptukey(abs(delta[i]) / se, k, nu, lower.tail = FALSE)), "\n")
}
# Scheffe: the critical value is the square root of (k - 1) times the F
# quantile, applied to the standard error of a difference between two means
crit <- sqrt((k - 1) * qf(0.95, k - 1, nu))
cat("Scheffe: Critical value =", six(crit), "\n")
sed <- sqrt(mse * (2 / n))
for (i in ord) {
tstat <- abs(delta[i]) / sed
cat(label[i], six(delta[i]), six(delta[i] - crit * sed), "to",
six(delta[i] + crit * sed), six(tstat),
pv(pf(tstat^2 / (k - 1), k - 1, nu, lower.tail = FALSE)), "\n")
}
# Newman-Keuls: the same statistic as Tukey's, but referred to the Studentized
# range for the number of ordered means that the two groups span (separation)
rnk <- rank(means)
for (i in ord) {
span <- abs(rnk[pairs[1, i]] - rnk[pairs[2, i]]) + 1
cat(label[i], six(delta[i]), span, six(abs(delta[i]) / se),
pv(ptukey(abs(delta[i]) / se, span, nu, lower.tail = FALSE)), "\n")
}
# Bonferroni: an ordinary t test of one planned comparison (Substance 1 with
# Substance 4) using the pooled variance, judged against a critical P of 1 minus
# the confidence level (0.05 here) divided by the number of comparisons; the
# simultaneous interval uses the same division. R's pairwise.t.test gives the P
# values so multiplied instead.
d14 <- means[1] - means[4]
tstat <- d14 / sed
comparisons <- 6
cat("Bonferroni: Mean A - Mean B =", six(d14), " Estimated std. error =", six(sed),
"\n")
cat("95% confidence interval =", six(d14 - qt(0.975, nu) * sed), "to",
six(d14 + qt(0.975, nu) * sed), "\n")
alpha <- 0.05 / comparisons
cat("Bonferroni-adjusted (simultaneous) 95% confidence interval for", comparisons,
"comparisons =", six(d14 - qt(1 - alpha / 2, nu) * sed), "to",
six(d14 + qt(1 - alpha / 2, nu) * sed),
paste0("(each comparison at ", formatC(100 * (1 - alpha), digits = 2, format = "f"),
"%)"), "\n")
cat("t =", six(tstat), " df =", nu, " ", pv(2 * pt(-abs(tstat), nu)), "\n")
cat("Bonferroni critical P for", comparisons, "comparisons =", six(alpha), "\n")
print(pairwise.t.test(y, group, p.adjust.method = "bonferroni"))
# Dunnett (each substance against Substance 1 as the control): base R has no
# Dunnett test (the multcomp package has one). Its critical value is the 95%
# point of the largest absolute value of the three t statistics, which share
# the control mean and so are correlated; with equal group sizes that
# probability is a double integral over the standard normal and the residual
# standard deviation, evaluated here numerically
m <- k - 1
chi <- function(u) exp(nu / 2 * log(nu) + (nu - 1) * log(u) - nu * u^2 / 2 -
(nu / 2 - 1) * log(2) - lgamma(nu / 2)) # density of s / sigma
inner <- function(u, d) sapply(u, function(ui) integrate(function(z) {
dnorm(z) * (pnorm(z + sqrt(2) * d * ui) - pnorm(z - sqrt(2) * d * ui))^m
}, -Inf, Inf, rel.tol = 1e-10)$value)
pmax_abs_t <- function(d) integrate(function(u) chi(u) * inner(u, d), 0, Inf,
rel.tol = 1e-9)$value
d <- uniroot(function(x) pmax_abs_t(x) - 0.95, c(1, 5), tol = 1e-9)$root
cat("Dunnett: Critical value (|d|) =", six(d), "\n")
for (j in c(4, 3, 2)) {
dj <- means[j] - means[1]
cat(names(means)[j], six(dj), six(dj - d * sed), "to", six(dj + d * sed),
pv(1 - pmax_abs_t(abs(dj) / sed)), "\n")
}