> 
> table(state,color)
      color
state  purple pink orange green
  ME        3    3      3     3
  VTNH      3    3      3     3
  MS        3    3      3     3
  WY        3    3      3     3
> 
> tapply(Y,student,mean)
        1         2         3         4         5         6         7         8 
 8.515436 15.239335  5.466240  4.620950 15.216931  9.449821 17.615383  4.918308 
        9        10        11        12 
12.258867 22.052906  5.233532  8.977321 
> tapply(Y,state,mean)
       ME      VTNH        MS        WY 
 3.609229  9.448168 12.513417 17.617528 
> tapply(Y,color,mean)
   purple      pink    orange     green 
 8.359209 12.754162 11.360606 10.714367 
> tapply(Y,square,mean)
        1         2         3 
 9.406910  9.932272 13.052075 
> 
> color.mod1 <- aov(Y~state+square/student+color)
> summary(color.mod1)
               Df  Sum Sq Mean Sq F value  Pr(>F)  
state           3 1235.39  411.80  4.1465 0.01428 *
square          2  124.25   62.12  0.6255 0.54180  
color           3  121.17   40.39  0.4067 0.74927  
square:student  9 1303.74  144.86  1.4586 0.20849  
Residuals      30 2979.34   99.31                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> TukeyHSD(color.mod1,"color")
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Y ~ state + square/student + color)

$color
                    diff        lwr       upr     p adj
pink-purple    4.3949527  -6.667476 15.457382 0.7040022
orange-purple  3.0013967  -8.061032 14.063825 0.8810919
green-purple   2.3551584  -8.707270 13.417587 0.9376464
orange-pink   -1.3935561 -12.455985  9.668873 0.9859020
green-pink    -2.0397943 -13.102223  9.022635 0.9581160
green-orange  -0.6462382 -11.708667 10.416191 0.9985401

>