> 
> mu <- c(244.3,321.7,243.3,286.0,300.7,251.0,249.7,299.7)
> s2 <- c(0.44,43.03,21.72,89.49,40.45,1.35,33.87,11.16)
> 
> y1 <- numeric(4)
> y  <- numeric(32)
> 
> for (i in 1:8) {
+ y1 <- rnorm(4,0,1)
+ y1 <- (y1-mean(y1))/sd(y1)
+ y[((i-1)*4+1):(i*4)] <- mu[i]+sqrt(s2[i])*y1
+ }
> 
> 
> a <- rep(c(-1,1,-1,1,1,-1,-1,1),each=4)
> b <- rep(c(-1,-1,1,1,-1,-1,1,1),each=4)
> c <- rep(c(-1,-1,-1,-1,1,1,1,1),each=4)
> 
> print(cbind(y,a,b,c))
             y  a  b  c
 [1,] 243.4201 -1 -1 -1
 [2,] 244.7590 -1 -1 -1
 [3,] 244.8616 -1 -1 -1
 [4,] 244.1593 -1 -1 -1
 [5,] 330.8149  1 -1 -1
 [6,] 319.9052  1 -1 -1
 [7,] 315.2119  1 -1 -1
 [8,] 320.8680  1 -1 -1
 [9,] 238.2248 -1  1 -1
[10,] 249.4653 -1  1 -1
[11,] 243.3866 -1  1 -1
[12,] 242.1234 -1  1 -1
[13,] 286.3971  1  1 -1
[14,] 293.4932  1  1 -1
[15,] 272.5408  1  1 -1
[16,] 291.5690  1  1 -1
[17,] 296.7262  1 -1  1
[18,] 304.0169  1 -1  1
[19,] 307.8965  1 -1  1
[20,] 294.1604  1 -1  1
[21,] 250.9584 -1 -1  1
[22,] 252.3079 -1 -1  1
[23,] 249.4905 -1 -1  1
[24,] 251.2432 -1 -1  1
[25,] 243.5448 -1  1  1
[26,] 255.3181 -1  1  1
[27,] 253.9696 -1  1  1
[28,] 245.9675 -1  1  1
[29,] 295.6806  1  1  1
[30,] 298.5380  1  1  1
[31,] 303.4200  1  1  1
[32,] 301.1614  1  1  1
> 
> partsize.mod1 <- aov(y~a*b*c)
> anova(partsize.mod1)
Analysis of Variance Table

Response: y
          Df  Sum Sq Mean Sq  F value    Pr(>F)    
a          1 24156.0 24156.0 800.1663 < 2.2e-16 ***
b          1   760.5   760.5  25.1915  3.96e-05 ***
c          1    16.8    16.8   0.5572 0.4626560    
a:b        1   591.7   591.7  19.5994 0.0001783 ***
a:c        1   208.1   208.1   6.8926 0.0148271 *  
b:c        1   591.7   591.7  19.5994 0.0001783 ***
a:b:c      1   612.5   612.5  20.2890 0.0001465 ***
Residuals 24   724.5    30.2                       
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> summary.lm(partsize.mod1)

Call:
aov(formula = y ~ a * b * c)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.4592  -2.2792   0.0225   3.4177   9.1149 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 274.5500     0.9713 282.666  < 2e-16 ***
a            27.4750     0.9713  28.287  < 2e-16 ***
b            -4.8750     0.9713  -5.019 3.96e-05 ***
c             0.7250     0.9713   0.746 0.462656    
a:b          -4.3000     0.9713  -4.427 0.000178 ***
a:c          -2.5500     0.9713  -2.625 0.014827 *  
b:c           4.3000     0.9713   4.427 0.000178 ***
a:b:c         4.3750     0.9713   4.504 0.000146 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 5.494 on 24 degrees of freedom
Multiple R-squared: 0.9738,     Adjusted R-squared: 0.9662 
F-statistic: 127.5 on 7 and 24 DF,  p-value: < 2.2e-16 

>