y <- c(6,4,5,3,4,3,6, 8,9,7,9,8,5,7, 6,7,7,7,8,5,7)
x <- c(3,1,3,1,2,1,4, 4,5,5,4,3,1,2, 3,2,2,3,4,1,4)
m <- gl(3,7)
plot(x[m==1], y[m==1], pch="A", xlim=c(0,6), ylim=c(2,10), xlab="Aptitude Scores", ylab="Achievement Scores")
points(x[m==2], y[m==2], pch="B")
points(x[m==3], y[m==3], pch="C")
anova(z <- lm(y~x+m))