> ## Chapter 4 > ## > data=read.table("table4=2.dat", skip=0, h=F) > data V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 1 - + + - - 7.78 7.78 7.81 7.7900 0.0003 -8.1117 2 + + + + - 8.15 8.18 7.88 8.0700 0.0273 -3.6009 3 - - + + - 7.50 7.56 7.50 7.5200 0.0012 -6.7254 4 + - + - - 7.59 7.56 7.75 7.6333 0.0104 -4.5627 5 - + - + - 7.94 8.00 7.88 7.9400 0.0036 -5.6268 6 + + - - - 7.69 8.09 8.06 7.9467 0.0496 -3.0031 7 - - - - - 7.56 7.62 7.44 7.5400 0.0084 -4.7795 8 + - - + - 7.56 7.81 7.69 7.6867 0.0156 -4.1583 9 - + + - + 7.50 7.25 7.12 7.2900 0.0373 -3.2888 10 + + + + + 7.88 7.88 7.44 7.7333 0.0645 -2.7406 11 - - + + + 7.50 7.56 7.50 7.5200 0.0012 -6.7254 12 + - + - + 7.63 7.75 7.56 7.6467 0.0092 -4.6849 13 - + - + + 7.32 7.44 7.44 7.4000 0.0048 -5.3391 14 + + - - + 7.56 7.69 7.62 7.6233 0.0042 -5.4648 15 - - - - + 7.18 7.18 7.25 7.2033 0.0016 -6.4171 16 + - - + + 7.81 7.50 7.59 7.6333 0.0254 -3.6717 > > # To avoid coding in R or Splus, 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 + } > > B=as.level2(data[,1]) > C=as.level2(data[,2]) > D=as.level2(data[,3]) > E=as.level2(data[,4]) > Q=as.level2(data[,5]) > > ybar = apply(data[,6:8], 1, mean) > lns2 = log( apply(data[,6:8], 1, var) ) > > # Table 4.5, note that factorial effects are twice the LS regression estimates, see sec 3.6 > g=lm(ybar~B+C+D+E+Q+B:Q+C:Q+D:Q+E:Q+B:C+B:D+C:D+B:C:Q+B:D:Q+B:E:Q) # BC=DE, BD=CE, BE=CD > round(2*g$coef,3) (Intercept) B C D E Q 15.272 0.221 0.176 0.029 0.104 -0.260 B:C B:D B:Q C:D C:Q D:Q 0.017 0.020 0.085 -0.035 -0.165 0.054 E:Q B:C:Q B:D:Q B:E:Q 0.027 0.010 -0.040 -0.047 > summary(g) Call: lm(formula = ybar ~ B + C + D + E + Q + B:Q + C:Q + D:Q + E:Q + B:C + B:D + C:D + B:C:Q + B:D:Q + B:E:Q) Residuals: ALL 16 residuals are 0: no residual degrees of freedom! Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.636042 B 0.110625 C 0.088125 D 0.014375 E 0.051875 Q -0.129792 B:C 0.008542 B:D 0.009792 B:Q 0.042292 C:D -0.017708 C:Q -0.082708 D:Q 0.026875 E:Q 0.013542 B:C:Q 0.005208 B:D:Q -0.020208 B:E:Q -0.023542 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 > > g2=lm(lns2~B+C+D+E+Q+B:Q+C:Q+D:Q+E:Q+B:C+B:D+C:D+B:C:Q+B:D:Q+B:E:Q) # BC=DE, BD=CE, BE=CD > round(2*g2$coef,3) (Intercept) B C D E Q -9.863 1.891 0.569 -0.247 0.216 0.280 B:C B:D B:Q C:D C:Q D:Q -0.002 0.425 -0.589 0.670 0.598 1.111 E:Q B:C:Q B:D:Q B:E:Q 0.129 -1.089 -0.432 0.854 > summary(g2) Call: lm(formula = lns2 ~ B + C + D + E + Q + B:Q + C:Q + D:Q + E:Q + B:C + B:D + C:D + B:C:Q + B:D:Q + B:E:Q) Residuals: ALL 16 residuals are 0: no residual degrees of freedom! Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -4.9313127 B 0.9454341 C 0.2843438 D -0.1237478 E 0.1077734 Q 0.1397582 B:C -0.0007902 B:D 0.2123448 B:Q -0.2943717 C:D 0.3352337 C:Q 0.2989006 D:Q 0.5553763 E:Q 0.0645709 B:C:Q -0.5446315 B:D:Q -0.2162354 B:E:Q 0.4267803 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 > > # we need the halfnormal function in halfnormal.s > source("halfnormal.s") > > # Fig 4.2 half-normal plot for location > halfnormal(2*g$coef[-1], label=T, type="n") # remove intercept > > # Fig 4.3 interaction plot > interaction.plot(Q, B, ybar, ylim=c(7.0, 8.0)) > interaction.plot(Q, C, ybar, ylim=c(7.0, 8.0)) > > # Fig 4.4 half-normal plot for dispersion > halfnormal(2*g2$coef[-1], label=T, type="n") # remove intercept > > # Fig 4.5 interaction plot > interaction.plot(Q, D, lns2) > > # Fig 4.6 interaction plot > interaction.plot(Q, 2*B+C, lns2) # quick plot, 3: B+C+, 1: B+C-, -1: B-C+, -3:B-C- > # another plot > B0=paste("B", data[,1], sep="") > C0=paste("C", data[,2], sep="") > B0C0=paste(B0,C0, sep="") > interaction.plot(Q, B0C0, lns2) >