> tapply(hairiness,list(twistlevel,testspeed),mean)
         1        2        3
1 744.6667 744.5000 775.8333
2 582.3333 584.8333 597.5000
3 492.3333 494.6667 507.1667


> tapply(hairiness,bobbin,mean)
       1        2        3        4        5        6 
632.5556 607.5556 608.5556 605.3333 614.4444 614.1111 


> 
> model1 <- aov(hairiness ~ twistlevel*testspeed + bobbin)
> summary(model1)
                     Df Sum Sq Mean Sq   F value    Pr(>F)    
twistlevel            2 611792  305896 1324.6920 < 2.2e-16 ***
testspeed             2   4637    2318   10.0402 0.0002928 ***
bobbin                5   4414     883    3.8231 0.0063403 ** 
twistlevel:testspeed  4    826     207    0.8946 0.4761743    
Residuals            40   9237     231                        
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


> 
> TukeyHSD(model1,"twistlevel")
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = hairiness ~ twistlevel * testspeed + bobbin)

$twistlevel
          diff       lwr        upr p adj
2-1 -166.77778 -179.1064 -154.44916     0
3-1 -256.94444 -269.2731 -244.61583     0
3-2  -90.16667 -102.4953  -77.83805     0



> TukeyHSD(model1,"testspeed")
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = hairiness ~ twistlevel * testspeed + bobbin)

$testspeed
         diff        lwr      upr     p adj
2-1  1.555556 -10.773060 13.88417 0.9494075
3-1 20.388889   8.060274 32.71750 0.0007072
3-2 18.833333   6.504718 31.16195 0.0017447