> 
> brand <- factor(brand)
> 
> rust.aov <- aov(Y~brand)
> anova(rust.aov)
Analysis of Variance Table

Response: Y
          Df Sum Sq Mean Sq F value    Pr(>F)    
brand      3  15954  5317.8  866.12 < 2.2e-16 ***
Residuals 36    221     6.1                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 
> TukeyHSD(rust.aov, "brand")
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Y ~ brand)

$brand
      diff        lwr         upr     p adj
2-1  46.30  43.315536  49.2844635 0.0000000
3-1  24.81  21.825536  27.7944635 0.0000000
4-1  -2.67  -5.654464   0.3144635 0.0933303
3-2 -21.49 -24.474464 -18.5055365 0.0000000
4-2 -48.97 -51.954464 -45.9855365 0.0000000
4-3 -27.48 -30.464464 -24.4955365 0.0000000

>