Score <- c(18.61,13.54,16.08,18.96,13.31,18.86,19.17,13.69,14.47,18.81,
18.22,19.42,20.25,25.25,20.36,22.43,17.22,22.31,19.58,23.96,
26.32,27.01,27.08,22.32,29.77)
Group <- as.factor(c(rep(1,5),rep(2,5),rep(3,5),rep(4,5),rep(5,5)))
Miller <- data.frame(Group, Score)Multiple Comparison Procedures - Miller Data
Creating the Data
This set of data is taken from Multiple Comparison Procedures by Larry E. Toothaker SAGE, Newbury Park, 1993, which he calls the Miller data as it is taken from Simultaneous Statistical Inference by R.G. Miller, Springer-Verlay, New York, 1981.
Creating the data in R:
library(dplyr)
Attaching package: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, uniongroup_by(Miller, Group) %>%
summarise(
count = n(),
mean = mean(Score, na.rm = TRUE),
sd = sd(Score, na.rm=TRUE)
)# A tibble: 5 × 4
Group count mean sd
<fct> <int> <dbl> <dbl>
1 1 5 16.1 2.68
2 2 5 17 2.68
3 3 5 20.7 2.68
4 4 5 21.1 2.68
5 5 5 26.5 2.68We can also visualise this data:
library(ggpubr)Loading required package: ggplot2ggboxplot(Miller, x="Group", y="Score",
color = "Group", palette= c("#4B0082","#5c3317", "#176608", "#990012", "#000099"),
order = c(1,2,3,4,5),
ylab = "Miller Score", xlab = "Treatment")
library(ggpubr)
ggline(Miller, x="Group", y="Score",
add = c("mean_se","jitter"),
order = c(1,2,3,4,5),
ylab = "Score", xlab = "Treatment")
Now you can carry out the one way ANOVA
results.anova <- aov(Score ~ Group, data = Miller)
summary(results.anova) Df Sum Sq Mean Sq F value Pr(>F)
Group 4 338.8 84.71 11.77 4.46e-05 ***
Residuals 20 144.0 7.20
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This shows that there is a significant result.
You can now carry our Tukey’s honest significant differences for the multiple pairwise comparisons.
TukeyHSD(results.anova) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Score ~ Group, data = Miller)
$Group
diff lwr upr p adj
2-1 0.9 -4.17823507 5.978235 0.9831073
3-1 4.6 -0.47823507 9.678235 0.0876742
4-1 5.0 -0.07823507 10.078235 0.0549250
5-1 10.4 5.32176493 15.478235 0.0000489
3-2 3.7 -1.37823507 8.778235 0.2272518
4-2 4.1 -0.97823507 9.178235 0.1517881
5-2 9.5 4.42176493 14.578235 0.0001553
4-3 0.4 -4.67823507 5.478235 0.9992543
5-3 5.8 0.72176493 10.878235 0.0203947
5-4 5.4 0.32176493 10.478235 0.0337298Alternatively you can carry out the corrected t-test using the Benjamini and Hochberg false discovery rate, Holm’s method or the Bonferroni correction.
pairwise.t.test(Miller$Score, Miller$Group,
p.adjust.method ="BH")
Pairwise comparisons using t tests with pooled SD
data: Miller$Score and Miller$Group
1 2 3 4
2 0.6686 - - -
3 0.0224 0.0517 - -
4 0.0160 0.0363 0.8161 -
5 5.5e-05 8.8e-05 0.0091 0.0117
P value adjustment method: BH pairwise.t.test(Miller$Score, Miller$Group,
p.adjust.method ="holm")
Pairwise comparisons using t tests with pooled SD
data: Miller$Score and Miller$Group
1 2 3 4
2 1.00000 - - -
3 0.06731 0.12407 - -
4 0.04790 0.10153 1.00000 -
5 5.5e-05 0.00016 0.02182 0.03279
P value adjustment method: holm pairwise.t.test(Miller$Score, Miller$Group,
p.adjust.method ="bonferroni")
Pairwise comparisons using t tests with pooled SD
data: Miller$Score and Miller$Group
1 2 3 4
2 1.00000 - - -
3 0.13463 0.41358 - -
4 0.07983 0.25381 1.00000 -
5 5.5e-05 0.00018 0.02728 0.04685
P value adjustment method: bonferroni In this case we could treat the first group as the control and compare the other groups to it. In this case we would use Dunnett’s test.
library(DescTools)
DunnettTest(x=Miller$Score, g=Miller$Group, control=5)
Dunnett's test for comparing several treatments with a control :
95% family-wise confidence level
$`5`
diff lwr.ci upr.ci pval
1-5 -10.4 -14.901421 -5.8985794 1.6e-05 ***
2-5 -9.5 -14.001421 -4.9985794 7.1e-05 ***
3-5 -5.8 -10.301421 -1.2985794 0.0095 **
4-5 -5.4 -9.901421 -0.8985794 0.0160 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Looking at the scatter-plot there is a good case to be made that you should carry out a test to see if there are differences between three sets of groups. They are 1 and 2 against 3 and 4 and against 5. To look at this I am going to use the linear model.
One-way ANOVA can be considered a linear model where you are testing to see if a single horizontal line passes through the means with an intercept but no gradient. If no such line passes through the means of all the groups then at least one group must have a different mean.
results.lm1 <- lm(Score ~ Group, data=Miller)
anova(results.lm1)Analysis of Variance Table
Response: Score
Df Sum Sq Mean Sq F value Pr(>F)
Group 4 338.84 84.71 11.765 4.463e-05 ***
Residuals 20 144.00 7.20
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(results.lm1)
Call:
lm(formula = Score ~ Group, data = Miller)
Residuals:
Min 1Q Median 3Q Max
-4.18 -2.48 -0.02 1.86 4.55
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 16.100 1.200 13.417 1.85e-11 ***
Group2 0.900 1.697 0.530 0.60172
Group3 4.600 1.697 2.711 0.01346 *
Group4 5.000 1.697 2.946 0.00798 **
Group5 10.400 1.697 6.128 5.47e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.683 on 20 degrees of freedom
Multiple R-squared: 0.7018, Adjusted R-squared: 0.6421
F-statistic: 11.77 on 4 and 20 DF, p-value: 4.463e-05This is the same ANOVA table as we had before. If we put group 5 as the reference then there is a significant difference between this group and all of the others in the linear model.
contrasts(Miller$Group) <- contr.treatment(5, base=5)
results.lm2 <- lm(Score ~ Group, data = Miller)
summary(results.lm2)
Call:
lm(formula = Score ~ Group, data = Miller)
Residuals:
Min 1Q Median 3Q Max
-4.18 -2.48 -0.02 1.86 4.55
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 26.500 1.200 22.083 1.62e-15 ***
Group1 -10.400 1.697 -6.128 5.47e-06 ***
Group2 -9.500 1.697 -5.598 1.77e-05 ***
Group3 -5.800 1.697 -3.418 0.00273 **
Group4 -5.400 1.697 -3.182 0.00468 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.683 on 20 degrees of freedom
Multiple R-squared: 0.7018, Adjusted R-squared: 0.6421
F-statistic: 11.77 on 4 and 20 DF, p-value: 4.463e-05I talked about a step-wise model where groups 1 and 2 and 3 and 4 are clustered together. You could calculate this using contrasts but combining groups is difficult and reading the results even harder. As a simpler alternative I just recoded the data with a dummy variable.
Dummy <- as.factor(c(rep(1,10),rep(3,10),rep(5,5)))
MillerD <- data.frame(Miller,Dummy)
results.anova1 <- aov(Score ~ Dummy, data = MillerD)
summary(results.anova1) Df Sum Sq Mean Sq F value Pr(>F)
Dummy 2 336.4 168.21 25.27 1.99e-06 ***
Residuals 22 146.4 6.66
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1contrasts(MillerD$Dummy) <- contr.treatment(3, base=3)
results.lm3 <- lm(Score ~ Dummy, data = MillerD)
summary(results.lm3)
Call:
lm(formula = Score ~ Dummy, data = MillerD)
Residuals:
Min 1Q Median 3Q Max
-4.18 -2.08 -0.18 2.26 4.35
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 26.500 1.154 22.969 < 2e-16 ***
Dummy1 -9.950 1.413 -7.042 4.59e-07 ***
Dummy2 -5.600 1.413 -3.963 0.00066 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.58 on 22 degrees of freedom
Multiple R-squared: 0.6967, Adjusted R-squared: 0.6692
F-statistic: 25.27 on 2 and 22 DF, p-value: 1.995e-06I can then compare this linear model to the Dunnett multiple comparison test.
DunnettTest(x=MillerD$Score, g=MillerD$Dummy, control="5")
Dunnett's test for comparing several treatments with a control :
95% family-wise confidence level
$`5`
diff lwr.ci upr.ci pval
1-5 -9.95 -13.245316 -6.654684 8.8e-07 ***
3-5 -5.60 -8.895316 -2.304684 0.0012 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1