> methgluc.mod1 <- aov(y ~ run + pressure)
> summary(methgluc.mod1)
            Df Sum Sq Mean Sq F value    Pr(>F)    
run          9 1394.7  154.96  5.0249  0.002529 ** 
pressure     4 3688.6  922.14 29.9020 3.026e-07 ***
Residuals   16  493.4   30.84                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> TukeyHSD(methgluc.mod1,"pressure")
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = y ~ run + pressure)

$pressure
             diff          lwr       upr     p adj
325-250 -2.444444 -12.26715022  7.378261 0.9377227
400-250  8.666667  -1.15603911 18.489372 0.0977009
475-250 15.277778   5.45507200 25.100484 0.0017011
550-250 25.166667  15.34396089 34.989372 0.0000062
400-325 11.111111   1.28840534 20.933817 0.0228901
475-325 17.722222   7.89951645 27.544928 0.0003836
550-325 27.611111  17.78840534 37.433817 0.0000019
475-400  6.611111  -3.21159466 16.433817 0.2827874
550-400 16.500000   6.67729423 26.322706 0.0008025
550-475  9.888889   0.06618311 19.711595 0.0480684