> # sec 3.11: Original Epitaxial Layer Growth Experiment : Table 3.8
> 
> data=read.table("originallayer.dat", skip=2, h=F)
> data
   V1 V2 V3 V4     V5     V6     V7     V8     V9    V10
1   -  -  -  + 14.812 14.774 14.772 14.794 14.860 14.914
2   -  -  -  - 13.768 13.778 13.870 13.896 13.932 13.914
3   -  -  +  + 14.722 14.736 14.774 14.778 14.682 14.850
4   -  -  +  - 13.860 13.876 13.932 13.846 13.896 13.870
5   -  +  -  + 14.886 14.810 14.868 14.876 14.958 14.932
6   -  +  -  - 14.182 14.172 14.126 14.274 14.154 14.082
7   -  +  +  + 14.758 14.784 15.054 15.058 14.938 14.936
8   -  +  +  - 13.996 13.988 14.044 14.028 14.108 14.060
9   +  -  -  + 15.272 14.656 14.258 14.718 15.198 15.490
10  +  -  -  - 14.324 14.092 13.536 13.588 13.964 14.328
11  +  -  +  + 13.918 14.044 14.926 14.962 14.504 14.136
12  +  -  +  - 13.614 13.202 13.704 14.264 14.432 14.228
13  +  +  -  + 14.648 14.350 14.682 15.034 15.384 15.170
14  +  +  -  - 13.970 14.448 14.326 13.970 13.738 13.738
15  +  +  +  + 14.184 14.402 15.544 15.424 15.036 14.470
16  +  +  +  - 13.866 14.130 14.256 14.000 13.640 13.592
> # 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[,1])
> B=as.level2(data[,2])
> C=as.level2(data[,3])
> D=as.level2(data[,4])
> 
> ybar = apply(data[,5:10], 1, mean)
> ybar
       1        2        3        4        5        6        7        8 
14.82100 13.85967 14.75700 13.88000 14.88833 14.16500 14.92133 14.03733 
       9       10       11       12       13       14       15       16 
14.93200 13.97200 14.41500 13.90733 14.87800 14.03167 14.84333 13.91400 
> s2 = apply(data[,5:10], 1, var)
> lns2 = log(s2)
> lns2
        1         2         3         4         5         6         7 
-5.770564 -5.311065 -5.703582 -6.984204 -5.916926 -5.484647 -4.106709 
        8         9        10        11        12        13        14 
-6.236717 -1.538078 -2.115847 -1.579159 -1.487046 -1.916322 -2.430302 
       15        16 
-1.118345 -2.653335 
> 
> # Table 3.9, note that factorial effects are twice the LS regression estimates, see sec 3.6
> g=lm(ybar~A*B*C*D)
> round(2*g$coef[-1],3)
      A       B       C       D     A:B     A:C     A:D     B:C     B:D 
 -0.055   0.142  -0.109   0.836  -0.032  -0.074  -0.025   0.047   0.010 
    C:D   A:B:C   A:B:D   A:C:D   B:C:D A:B:C:D 
 -0.037   0.060   0.067  -0.056   0.098   0.036 
> 
> g2=lm(lns2~A*B*C*D)
> round(2*g2$coef[-1],3)
      A       B       C       D     A:B     A:C     A:D     B:C     B:D 
  3.834   0.078   0.077   0.632  -0.428   0.214   0.002   0.331   0.305 
    C:D   A:B:C   A:B:D   A:C:D   B:C:D A:B:C:D 
  0.582  -0.335   0.086  -0.494   0.314   0.109 
> 
> # we create a function for convenience
> 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))
+ }
> 
> # Fig 3.11 half-normal plot for location
> halfnormal(2*g$coef[-1]) # remove intercept
> halfnormal(2*g$coef[-1], label=T, xlim=c(-.1, 2.5))  # add labels
> # Fig 3.11 half-normal plot for dispersion 
> halfnormal(2*g2$coef[-1]) # remove intercept
> halfnormal(2*g2$coef[-1], label=T, xlim=c(-.1, 2.5))  # add labels