> # Sec 3.1-3.10: Adapted Epitaxial Layer Growth Experiment
> # Table 3.1, 
> data=read.table("adaptedlayer.dat", skip=2, h=F)
> data
   V1 V2 V3 V4 V5     V6     V7     V8     V9    V10    V11
1   1  -  -  -  + 14.506 14.153 14.134 14.339 14.953 15.455
2   2  -  -  -  - 12.886 12.963 13.669 13.869 14.145 14.007
3   3  -  -  +  + 13.926 14.052 14.392 14.428 13.568 15.074
4   4  -  -  +  - 13.758 13.992 14.808 13.554 14.283 13.904
5   5  -  +  -  + 14.629 13.940 14.466 14.538 15.281 15.046
6   6  -  +  -  - 14.059 13.989 13.666 14.706 13.863 13.357
7   7  -  +  +  + 13.800 13.896 14.887 14.902 14.461 14.454
8   8  -  +  +  - 13.707 13.623 14.210 14.042 14.881 14.378
9   9  +  -  -  + 15.050 14.361 13.916 14.431 14.968 15.294
10 10  +  -  -  - 14.249 13.900 13.065 13.143 13.708 14.255
11 11  +  -  +  + 13.327 13.457 14.368 14.405 13.932 13.552
12 12  +  -  +  - 13.605 13.190 13.695 14.259 14.428 14.223
13 13  +  +  -  + 14.274 13.904 14.317 14.754 15.188 14.923
14 14  +  +  -  - 13.775 14.586 14.379 13.775 13.382 13.382
15 15  +  +  +  + 13.723 13.914 14.913 14.808 14.469 13.973
16 16  +  +  +  - 14.031 14.467 14.675 14.252 13.658 13.578

> # Since PC and Mac have different coding for "-" and "+", we create a function to code -/+ as -1/+1
> as.level2=function(x)
+ {
+ 	y = rep(0, length(x))
+ 	y[ x=="+" ] = 1
+ 	y[ x=="-" ] = -1
+ 	y
+ }
> 
> A=as.level2(data[,2])
> B=as.level2(data[,3])
> C=as.level2(data[,4])
> D=as.level2(data[,5])
> A
 [1] -1 -1 -1 -1 -1 -1 -1 -1  1  1  1  1  1  1  1  1
> 
 
> ybar = apply(data[,6:11], 1, mean)
> ybar
       1        2        3        4        5        6        7        8 
14.59000 13.58983 14.24000 14.04983 14.65000 13.94000 14.40000 14.14017 
       9       10       11       12       13       14       15       16 
14.67000 13.72000 13.84017 13.90000 14.56000 13.87983 14.30000 14.11017 
> s2 = apply(data[,6:11], 1, var)
> lns2 = log(s2)
> lns2
        1         2         3         4         5         6         7 
-1.310107 -1.234610 -1.317121 -1.624742 -1.508762 -1.585514 -1.505235 
        8         9        10        11        12        13        14 
-1.537111 -1.314012 -1.301539 -1.514787 -1.473908 -1.482832 -1.373975 
       15        16 
-1.385296 -1.649836 
> 
> # Fig 3.5
> par(mfrow=c(2,3))
> interaction.plot(B, A, ybar, ylim=c(13.6, 14.6))
> interaction.plot(C, A, ybar, ylim=c(13.6, 14.6))
> interaction.plot(D, A, ybar, ylim=c(13.6, 14.6))
> interaction.plot(C, B, ybar, ylim=c(13.6, 14.6))
> interaction.plot(D, B, ybar, ylim=c(13.6, 14.6))
> interaction.plot(D, C, ybar, ylim=c(13.6, 14.6))
> 
> # Fig 3.6
> interaction.plot(D, C, ybar, ylim=c(13.6, 14.6))
> interaction.plot(C, D, ybar, ylim=c(13.6, 14.6))
> 
> par(mfrow=c(1,1))  # reset output device
> 
> 
> # Table 3.6
> g=lm(ybar~A*B*C*D)
> summary(g)

Call:
lm(formula = ybar ~ A * B * C * D)

Residuals:
ALL 16 residuals are 0: no residual degrees of freedom!

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.161250                            
A           -0.038729                            
B            0.086271                            
C           -0.038708                            
D            0.245021                            
A:B          0.003708                            
A:C         -0.046229                            
A:D         -0.025000                            
B:C          0.028771                            
B:D         -0.015042                            
C:D         -0.172521                            
A:B:C        0.048750                            
A:B:D        0.012521                            
A:C:D       -0.015000                            
B:C:D        0.054958                            
A:B:C:D      0.009979                            

Residual standard error: NaN on 0 degrees of freedom
Multiple R-Squared:     1,	Adjusted R-squared:   NaN 
F-statistic:   NaN on 15 and 0 DF,  p-value:   NaN 

> # factorial effects are twice the LS regression estimates, see sec 3.6
> 2*g$coef
 (Intercept)            A            B            C            D          A:B 
28.322500000 -0.077458333  0.172541667 -0.077416667  0.490041667  0.007416667 
         A:C          A:D          B:C          B:D          C:D        A:B:C 
-0.092458333 -0.050000000  0.057541667 -0.030083333 -0.345041667  0.097500000 
       A:B:D        A:C:D        B:C:D      A:B:C:D 
 0.025041667 -0.030000000  0.109916667  0.019958333 
> 
> g2=lm(lns2~A*B*C*D)
> 2*g2$coef
  (Intercept)             A             B             C             D 
-2.8899232671  0.0158770768 -0.1172170200 -0.1120855148  0.0553852385 
          A:B           A:C           A:D           B:C           B:D 
 0.0452936809 -0.0257819863 -0.0298029198  0.0804866861  0.0106925625 
          C:D         A:B:C         A:B:D         A:C:D         B:C:D 
 0.0854040848 -0.0317823268  0.0415664954  0.0008436362 -0.0032738955 
      A:B:C:D 
 0.1037245505 
> 
> # Fig 3.8 normal plot, 
> # we create a function for convenience
> normalplot=function(y, label=F, ...)
+ {
+ 	n=length(y)
+ 	x=seq(0.5/n, 1.0-0.5/n, by=1/n)
+ 	x = qnorm(x)
+ 	y = sort(y)
+ 	qqplot(x, y, xlab="normal quantiles", ylab="effects", ...)
+ 	if(label) text(x,y, names(y))
+ }
> 
> normalplot(2*g$coef[-1] ) # remove intercept
> normalplot(2*g$coef[-1], label=T, xlim=c(-2.1, 2.1) ) # add labels
> # or equivalently
> qqnorm(2*g$coef[-1]) # remove intercept
> 
> # Fig 3.9 half-normal plot, without labeling
> halfnormal=function(y, label=F, ...)
+ {
+ 	n=length(y)
+ 	x=seq(0.5+0.25/n, 1.0-0.25/n, by=0.5/n)
+ 	x = qnorm(x)
+ 	y = sort(abs(y))
+ 	qqplot(x, y, xlab="half-normal quantiles", ylab="absolute effects", ...)
+ 	if(label) text(x,y, names(y))
+ }
> 
> halfnormal(2*g2$coef[-1]) # remove intercept
> halfnormal(2*g2$coef[-1], label=T, xlim=c(-.1, 2.5))  # add labels
>