Enter data

The following example comes from the book Statistical Models in S (Chambers and Hastie, 1992). Unfortunately I had to enter the data by hand because I couldn’t find an electronic version.

Method <- factor(rep(c("I", "I", "II", "II"), 9))
Physique <- factor(rep(c("slight", "average", "heavy"), each = 12),
                   levels = c("slight", "average", "heavy"))
Team <- factor(rep(1:9, each = 4))
guns <- data.frame(Method, Physique, Team)
guns$Rounds <- c(20.2, 24.1, 14.2, 16.2,
                 26.2, 26.9, 18.0, 19.1,
                 23.8, 24.9, 12.5, 15.4,
                 22.0, 23.5, 14.1, 16.1,
                 22.6, 24.6, 14.0, 18.1,
                 22.9, 25.0, 13.7, 16.0,
                 23.1, 22.9, 14.1, 16.1, 
                 22.9, 23.7, 12.2, 13.8, 
                 21.8, 23.5, 12.7, 15.1)

About the data

The data come from another book called Fundamental Concepts in the Design of Experiments (Hicks and Turner, 1991). It’s based on an experiment about firing naval guns. Two different reloading methods were tested by gunners corresponding to three different Physiques (slight, average, heavy). Three teams were randomly chosen to represent each of the three physique groupings. Each team was then tested twice using each method for a total of 36 trials. The response was the number of rounds fired per minute.

summary(guns)
##  Method     Physique       Team        Rounds     
##  I :18   slight :12   1      : 4   Min.   :12.20  
##  II:18   average:12   2      : 4   1st Qu.:14.88  
##          heavy  :12   3      : 4   Median :19.65  
##                       4      : 4   Mean   :19.33  
##                       5      : 4   3rd Qu.:23.50  
##                       6      : 4   Max.   :26.90  
##                       (Other):12

Notice how Team is nested within Physique in the plot below. The Team variable is only meaningful within each Physique.

library(ggplot2)
ggplot(guns) +
  aes(x = Team, y = Rounds, color = Method) +
  geom_point() +
  facet_wrap(~ Physique, labeller = "label_both", scales = "free_x") +
  labs(title = "Team is nested within Physique")

Therefore we nest Team within Physique when we perform the ANOVA:

m <- aov(Rounds ~ Method + Physique/Team, guns)
summary(m)
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## Method         1  652.0   652.0  316.84 4.41e-16 ***
## Physique       2   16.1     8.0    3.90   0.0330 *  
## Physique:Team  6   39.3     6.5    3.18   0.0178 *  
## Residuals     26   53.5     2.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Most of the variability in Rounds is due to Method. Notice the huge Sum Sq value. Method I seems to be far superior to Method II. We might follow up the ANOVA with a comparison between Methods. It appears Method I produces about 8 more rounds on average.

TukeyHSD(m, which = "Method")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Rounds ~ Method + Physique/Team, data = guns)
## 
## $Method
##           diff       lwr       upr p adj
## II-I -8.511111 -9.493963 -7.528259     0

We might also want to compare Physique. There appears to be a significant difference between Heavy and Slight physiques, where slight performed better. But that’s probably due to Team 2 who seemed to perform exceptionally well.

TukeyHSD(m, which = "Physique")
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Rounds ~ Method + Physique/Team, data = guns)
## 
## $Physique
##                      diff       lwr        upr     p adj
## average-slight -0.7416667 -2.196851  0.7135174 0.4262335
## heavy-slight   -1.6333333 -3.088517 -0.1781493 0.0255522
## heavy-average  -0.8916667 -2.346851  0.5635174 0.2970760