pdf("hairyarn.pdf")

hairyarn <- read.fwf("C:\\data\\hairyarn.txt", width=c(8,8,8,8), 
col.names=c("twistlevel","testspeed","bobbin","hairiness"))

attach(hairyarn)

twistlevel <- factor(twistlevel)
testspeed <- factor(testspeed)
bobbin <- factor(bobbin)

tapply(hairiness,list(twistlevel,testspeed),mean)
tapply(hairiness,bobbin,mean)

model1 <- aov(hairiness ~ twistlevel*testspeed + bobbin)
summary(model1)

TukeyHSD(model1,"twistlevel")
TukeyHSD(model1,"testspeed")


dev.off()