> ## Layer Growth Experiment: sec. 10.5.1 > # > data=read.table("http://www.stat.ucla.edu/~hqxu/stat225/data/layerrobust.dat") > data V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 1 -1 -1 -1 1 -1 -1 -1 -1 15.3182 15.4279 15.2657 15.4056 14.2908 14.1924 14.2714 14.1876 2 -1 -1 -1 1 1 1 1 1 14.9306 14.8954 14.9210 15.1349 14.8030 14.7193 14.6960 14.7635 3 -1 -1 1 -1 -1 -1 1 1 14.0121 13.9386 14.2118 14.0789 13.8793 13.9213 13.8532 14.0849 4 -1 -1 1 -1 1 1 -1 -1 14.2444 14.2573 14.3951 14.3724 13.4054 13.4788 13.5878 13.5167 5 -1 1 -1 -1 -1 1 -1 1 14.1492 14.1654 14.1487 14.2765 14.1736 14.0360 14.1398 14.0796 6 -1 1 -1 -1 1 -1 1 -1 14.2204 14.3028 14.2689 14.4104 13.2539 13.3338 13.1920 13.4430 7 -1 1 1 1 -1 1 1 -1 15.2969 15.5209 15.4200 15.2077 14.0623 14.0888 14.1766 14.0528 8 -1 1 1 1 1 -1 -1 1 15.0100 15.0618 15.5724 15.4668 14.3068 14.4055 14.6780 14.5811 9 1 -1 -1 -1 -1 1 1 -1 14.9039 14.7952 14.1886 14.6254 13.7259 13.2934 12.6502 13.2666 10 1 -1 -1 -1 1 -1 -1 1 13.7546 14.3229 14.2224 13.8209 13.8953 14.5597 14.4492 13.7064 11 1 -1 1 1 -1 1 -1 1 14.1936 14.4295 15.5537 15.2200 14.2201 14.3974 15.2757 15.0363 12 1 -1 1 1 1 -1 1 -1 14.5640 14.4670 15.2293 15.1099 13.5228 13.5828 14.2822 13.8449 13 1 1 -1 1 -1 -1 1 1 14.7437 14.1827 14.9695 15.5484 14.5335 14.2492 14.6710 15.2799 14 1 1 -1 1 1 1 -1 -1 15.8717 15.2239 14.9700 16.0001 14.5676 14.0310 13.7099 14.6375 15 1 1 1 -1 -1 -1 -1 -1 14.2537 13.8368 14.1332 15.1681 12.9012 12.7071 13.1484 13.8940 16 1 1 1 -1 1 1 1 1 13.8136 14.0745 14.4313 13.6862 13.9532 14.0830 14.1119 13.5963 > # change some variable names > names(data)[1:8] = LETTERS[1:8] > names(data) [1] "A" "B" "C" "D" "E" "F" "G" "H" "V9" "V10" "V11" "V12" "V13" [14] "V14" "V15" "V16" > > # responses for location and dispersion models > ybar = apply(data[,9:16], 1, mean) > lns2 = log( apply(data[,9:16], 1, var) ) > > # make A-H be available > attach(data) > > # we need the halfnormal function in halfnormal.R > source("http://www.stat.ucla.edu/~hqxu/stat225/R/halfnormal.R") > > # we need to know the aliasing among factorial effects > # D=-ABC, F=ABE, G=ACE, H=BCE > # we can entertain 8 ME + 7 2fi's invloving A > > # Fig 10.2 half-normal plot for location > g=lm(ybar~A+B+C+D+E+F+G+H+A:B+A:C+A:D+A:E+A:F+A:G+A:H) > g # summary(g) Call: lm(formula = ybar ~ A + B + C + D + E + F + G + H + A:B + A:C + A:D + A:E + A:F + A:G + A:H) Coefficients: (Intercept) A B C D E 14.351997 -0.038061 0.014948 -0.057183 0.401917 -0.012541 F G H A:B A:C A:D 0.048841 -0.053950 0.086762 0.014306 -0.046700 -0.024641 A:E A:F A:G A:H 0.014105 0.029011 -0.010611 0.005102 > halfnormal(2*g$coef[-1], label=T, type="n") # remove intercept > > # Fig 10.2 half-normal plot for dispersion > g=lm(lns2~A+B+C+D+E+F+G+H+A:B+A:C+A:D+A:E+A:F+A:G+A:H) > g # summary(g) Call: lm(formula = lns2 ~ A + B + C + D + E + F + G + H + A:B + A:C + A:D + A:E + A:F + A:G + A:H) Coefficients: (Intercept) A B C D E -1.82210 0.61910 0.10219 0.16566 0.42621 0.02931 F G H A:B A:C A:D -0.20828 -0.10912 -0.98162 -0.13783 -0.25271 -0.22503 A:E A:F A:G A:H -0.35173 0.24258 -0.03052 0.30027 > halfnormal(2*g$coef[-1], label=T, type="n") # remove intercept > > > # sec. 10.5.2 Response modeling > # > # we have to create L, Ml, Mq, Mc, and modify the layout of the design > # > # arrange the response in a column-wise fashion > > data=read.table("layerrobust.dat") > data[,4] = - data[,4] # column D has a wrong sign. change it > attach(data) > > y = c(V9, V10, V11, V12, V13, V14, V15, V16) > > # The book codes bottom as +1 and top as -1, which is inconsistent with the level settings in Table 10.1. > L = rep(c(1, -1), c(64,64)) > L [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 [26] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 [51] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 [76] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 [101] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 [126] -1 -1 -1 > M = rep( rep(1:4, c(16,16,16,16)), 2) > M [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 [39] 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 1 1 1 1 1 1 1 1 1 1 1 1 [77] 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 [115] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 > > Ml = rep( rep(c(1,1,-1,-1), c(16,16,16,16)), 2) > Mq = rep( rep(c(1,-1,-1,1), c(16,16,16,16)), 2) > Mc = rep( rep(c(1,-1,1,-1), c(16,16,16,16)), 2) > > # repeat A-H 8 times > A = rep(V1, 8) > B = rep(V2, 8) > C = rep(V3, 8) > D = rep(V4, 8) > E = rep(V5, 8) > F = rep(V6, 8) > G = rep(V7, 8) > H = rep(V8, 8) > > ## There are 16*8=128 observations, so we can entertain 127 factorial effects. > g = lm(y ~ (A+B+C+D+E+F+G+H+A:B+A:C+A:D+A:E+A:F+A:G+A:H) * (L*(Ml+Mq+Mc)) ) > halfnormal(2*g$coef[-1], label=T, type="n") # remove intercept > > ## equation (10.10) > g #summary(g) Call: lm(formula = y ~ (A + B + C + D + E + F + G + H + A:B + A:C + A:D + A:E + A:F + A:G + A:H) * (L * (Ml + Mq + Mc))) Coefficients: (Intercept) A B C D E 1.435e+01 -3.806e-02 1.495e-02 -5.718e-02 4.019e-01 -1.254e-02 F G H L Ml Mq 4.884e-02 -5.395e-02 8.676e-02 3.295e-01 -9.016e-02 2.736e-02 Mc A:B A:C A:D A:E A:F -1.496e-02 1.431e-02 -4.670e-02 -2.464e-02 1.410e-02 2.901e-02 A:G A:H L:Ml L:Mq L:Mc A:L -1.061e-02 5.102e-03 -2.313e-02 4.028e-03 -4.609e-03 -8.820e-03 A:Ml A:Mq A:Mc B:L B:Ml B:Mq -4.707e-02 2.931e-02 -4.758e-03 3.498e-02 -1.383e-02 5.694e-02 B:Mc C:L C:Ml C:Mq C:Mc D:L -1.283e-02 7.902e-03 -8.307e-02 -5.478e-02 2.327e-02 4.572e-02 D:Ml D:Mq D:Mc E:L E:Ml E:Mq -6.653e-02 2.225e-02 7.861e-03 -1.188e-02 2.936e-02 -4.780e-02 E:Mc F:L F:Ml F:Mq F:Mc G:L 1.409e-02 -1.681e-03 5.823e-02 2.666e-03 1.073e-02 7.000e-04 G:Ml G:Mq G:Mc H:L H:Ml H:Mq 3.173e-02 -1.878e-03 1.356e-02 -2.389e-01 -4.190e-02 -6.822e-02 H:Mc A:L:Ml A:L:Mq A:L:Mc B:L:Ml B:L:Mq 2.238e-02 -9.817e-03 9.533e-03 1.995e-03 3.677e-03 4.453e-04 B:L:Mc C:L:Ml C:L:Mq C:L:Mc D:L:Ml D:L:Mq 5.533e-03 4.844e-05 -2.692e-03 1.001e-02 -3.217e-03 6.098e-03 D:L:Mc E:L:Ml E:L:Mq E:L:Mc F:L:Ml F:L:Mq -3.045e-03 -9.725e-03 1.001e-02 -3.781e-04 -7.344e-04 -1.728e-02 F:L:Mc G:L:Ml G:L:Mq G:L:Mc H:L:Ml H:L:Mq -5.312e-03 -5.316e-03 -1.483e-02 1.060e-02 -1.935e-02 -2.594e-03 H:L:Mc A:B:L A:B:Ml A:B:Mq A:B:Mc A:C:L 7.650e-03 -1.716e-02 -2.978e-03 7.125e-02 -1.168e-02 -2.838e-02 A:C:Ml A:C:Mq A:C:Mc A:D:L A:D:Ml A:D:Mq -5.223e-02 -2.582e-02 1.064e-02 -4.032e-02 -6.367e-02 3.575e-02 A:D:Mc A:E:L A:E:Ml A:E:Mq A:E:Mc A:F:L 5.406e-04 -2.667e-02 6.035e-02 -5.280e-02 2.444e-02 3.802e-02 A:F:Ml A:F:Mq A:F:Mc A:G:L A:G:Ml A:G:Mq 3.581e-02 6.255e-03 1.011e-03 1.255e-02 1.165e-02 -9.277e-03 A:G:Mc A:H:L A:H:Ml A:H:Mq A:H:Mc A:B:L:Ml 2.520e-02 -5.218e-02 -1.376e-02 -8.163e-02 1.605e-02 1.672e-03 A:B:L:Mq A:B:L:Mc A:C:L:Ml A:C:L:Mq A:C:L:Mc A:D:L:Ml 5.116e-03 7.934e-03 -7.563e-04 3.928e-03 -1.031e-02 -4.206e-03 A:D:L:Mq A:D:L:Mc A:E:L:Ml A:E:L:Mq A:E:L:Mc A:F:L:Ml 4.791e-03 1.160e-02 4.077e-03 3.484e-04 -6.352e-03 -8.714e-03 A:F:L:Mq A:F:L:Mc A:G:L:Ml A:G:L:Mq A:G:L:Mc A:H:L:Ml -2.138e-02 7.883e-03 -8.139e-03 -3.177e-03 -4.858e-03 -8.922e-04 A:H:L:Mq A:H:L:Mc -5.636e-03 1.077e-03 > > # the most sign't 7 factorial effects > round( sort(abs(g$coef)), 3)[121:128] A:H:Mq C:Ml H Ml H:L L 0.082 0.083 0.087 0.090 0.239 0.330 D (Intercept) 0.402 14.352 > > # Fig 10.5-10.7 > par(mfrow=c(2,2)) > interaction.plot(L, H, y) > interaction.plot(M, C, y) > > interaction.plot(M, 2*A+H, y) # quick plot, 3: A+H+, 1: A+H-, -1: A-H+, -3:A-H- > interaction.plot(M, D, y) >