Asphalt Data from the textbook Example 11.7

Prepare the dataset.

data = data.frame(conductivity = c(.835, .845, .822, .826, .785, .795, 
                                   .855, .865, .832, .836, .790, .800, 
                                   .815, .825, .800, .820, .770, .790 ),
                  
                  AsphGr = as.factor(rep(c("PG58", "PG64", "PG70"), each=6)),
                  AggCont = as.factor(rep(rep(c("38", "41", "44"), each=2), 3))
                  )

head(data)

Run the two-way ANOVA without interaction term.

fit1 = aov(conductivity~AsphGr + AggCont, data=data)
print(summary(fit1))
##             Df   Sum Sq  Mean Sq F value   Pr(>F)    
## AsphGr       2 0.002089 0.001045    13.7  0.00063 ***
## AggCont      2 0.008297 0.004149    54.4 4.83e-07 ***
## Residuals   13 0.000991 0.000076                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We can also add the interaction term.

fit2 = aov(conductivity~AsphGr + AggCont + AsphGr*AggCont, data=data)
print(summary(fit2))
##                Df   Sum Sq  Mean Sq F value   Pr(>F)    
## AsphGr          2 0.002089 0.001045  14.117  0.00168 ** 
## AggCont         2 0.008297 0.004149  56.063 8.31e-06 ***
## AsphGr:AggCont  4 0.000325 0.000081   1.099  0.41356    
## Residuals       9 0.000666 0.000074                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Remark: the interaction term can also be AsphGr:AggCont.

Built-in Dataset ToothGrowth

Check the first few lines in the dataset

head(ToothGrowth)

Now we can run a two-way ANOVA of len against supp and dose along with the interaction.

fit3 = aov(len~supp + dose + supp:dose, data=ToothGrowth)
print(summary(fit3))
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## supp         1  205.4   205.4  12.317 0.000894 ***
## dose         1 2224.3  2224.3 133.415  < 2e-16 ***
## supp:dose    1   88.9    88.9   5.333 0.024631 *  
## Residuals   56  933.6    16.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Three-way ANOVA

We use the CO2 dataset with a bit modification on the Plant

data = CO2
data$Plant = as.factor(rep(rep(c("1", "2", "3"), each=7), 4))
head(data)

It has three factors: Plant, Type and Treatment. The outcome is uptake. We consider a three-way ANOVA without interactions.

fit4 = aov(uptake ~ ., data=data)
print(summary(fit4))
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Plant        2     20      10   0.257    0.774    
## Type         1   3366    3366  86.119 3.09e-14 ***
## Treatment    1    988     988  25.284 3.08e-06 ***
## conc         1   2285    2285  58.469 4.54e-11 ***
## Residuals   78   3048      39                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Add the two-way interaction.

fit5 = aov(uptake ~ .^2, data=data)
print(summary(fit5))
##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## Plant            2     20      10   0.285   0.7526    
## Type             1   3366    3366  95.654 1.16e-14 ***
## Treatment        1    988     988  28.084 1.32e-06 ***
## conc             1   2285    2285  64.943 1.56e-11 ***
## Plant:Type       2    112      56   1.587   0.2118    
## Plant:Treatment  2     41      20   0.579   0.5632    
## Plant:conc       2      2       1   0.035   0.9655    
## Type:Treatment   1    226     226   6.416   0.0136 *  
## Type:conc        1    208     208   5.912   0.0176 *  
## Treatment:conc   1     32      32   0.906   0.3445    
## Residuals       69   2428      35                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Textbook Example 11.9

data = data.frame(viberation = c(13.1, 13.2, 16.3, 15.8, 13.7, 14.3, 15.7, 15.8, 13.5, 12.5, 
                                 15.0, 14.8, 15.7, 16.4, 13.9, 14.3, 13.7, 14.2, 13.4, 13.8, 
                                 14.0, 14.3, 17.2, 16.7, 12.4, 12.3, 14.4, 13.9, 13.2, 13.1),
                  
                  casmater = as.factor(rep(c("Steel", "Aluminum", "Plastic"), each=10)),
                  source = as.factor(rep(rep(1:5, each=2), 3))
                  )
print(data)
##    viberation casmater source
## 1        13.1    Steel      1
## 2        13.2    Steel      1
## 3        16.3    Steel      2
## 4        15.8    Steel      2
## 5        13.7    Steel      3
## 6        14.3    Steel      3
## 7        15.7    Steel      4
## 8        15.8    Steel      4
## 9        13.5    Steel      5
## 10       12.5    Steel      5
## 11       15.0 Aluminum      1
## 12       14.8 Aluminum      1
## 13       15.7 Aluminum      2
## 14       16.4 Aluminum      2
## 15       13.9 Aluminum      3
## 16       14.3 Aluminum      3
## 17       13.7 Aluminum      4
## 18       14.2 Aluminum      4
## 19       13.4 Aluminum      5
## 20       13.8 Aluminum      5
## 21       14.0  Plastic      1
## 22       14.3  Plastic      1
## 23       17.2  Plastic      2
## 24       16.7  Plastic      2
## 25       12.4  Plastic      3
## 26       12.3  Plastic      3
## 27       14.4  Plastic      4
## 28       13.9  Plastic      4
## 29       13.2  Plastic      5
## 30       13.1  Plastic      5

Then we run a mixed-effect ANOVA by treating the supply as a random effect.

fit6 = aov(viberation ~ casmater*source + Error(source), data=data)
print(summary(fit6))
## 
## Error: source
##        Df Sum Sq Mean Sq
## source  4  36.67   9.169
## 
## Error: Within
##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## casmater         2  0.705  0.3523   3.165   0.0713 .  
## casmater:source  8 11.605  1.4507  13.030 1.76e-05 ***
## Residuals       15  1.670  0.1113                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1