> dat0=read.table("http://www2.isye.gatech.edu/%7Ejeffwu/book/data/seatbelt.dat", h=F) > dimnames(dat0)[[2]]=c("run", "A", "B", "C", "D", "strength1", "strength2", "strength3", "flash1", "flash2", "flash3") > dat=dat0 > dat3=cbind(rbind(dat[,1:5],dat[,1:5],dat[,1:5]), strength=c(dat[,6], dat[,7], dat[,8]), flash=c(dat[,9], dat[,10], dat[,11])) > source("http://www.stat.ucla.edu/~hqxu/stat225/R/halfnormal.R") > library(MASS) > > # location model > dat=dat3 > y=dat3$strength # strength > # y=dat3$flash # flash > > # Note: D=A+B+C (mod 3) > dat[,5] - (dat[,2]+dat[,3]+dat[,4]) %%3 [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [53] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 > > A=as.ordered(dat[,2]) > B=as.ordered(dat[,3]) > C=as.ordered(dat[,4]) > D=as.ordered(dat[,5]) > > > AB=as.ordered((dat[,2]+dat[,3]) %%3) > AB2=as.ordered((dat[,2]+2*dat[,3]) %%3) > AC=as.ordered((dat[,2]+dat[,4]) %%3) > AC2=as.ordered((dat[,2]+2*dat[,4]) %%3) > BC=as.ordered((dat[,3]+dat[,4]) %%3) > BC2=as.ordered((dat[,3]+2*dat[,4]) %%3) > ABC=as.ordered((dat[,2]+dat[,3]+dat[,4]) %%3) > ABC2=as.ordered((dat[,2]+dat[,3]+2*dat[,4]) %%3) > AB2C=as.ordered((dat[,2]+2*dat[,3]+dat[,4]) %%3) > AB2C2=as.ordered((dat[,2]+2*dat[,3]+2*dat[,4]) %%3) > > # AB=CD2, AC=BD2, BC=AD2 > AD=as.ordered((dat[,2]+dat[,5]) %%3) > AD2=as.ordered((dat[,2]+2*dat[,5]) %%3) > BD=as.ordered((dat[,3]+dat[,5]) %%3) > BD2=as.ordered((dat[,3]+2*dat[,5]) %%3) > CD=as.ordered((dat[,4]+dat[,5]) %%3) > CD2=as.ordered((dat[,4]+2*dat[,5]) %%3) # =AB > > > g=lm(y~(A+B+C)^3) > anova(g) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) A 2 34621746 17310873 85.5828 < 2.2e-16 *** B 2 938539 469270 2.3200 0.107992 C 2 9549481 4774741 23.6057 4.3e-08 *** A:B 4 3298246 824561 4.0765 0.005846 ** A:C 4 3872179 968045 4.7859 0.002231 ** B:C 4 448348 112087 0.5541 0.696829 A:B:C 8 5206919 650865 3.2178 0.004620 ** Residuals 54 10922599 202270 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > g=lm(y~A+B+AB+AB2+C+AC+AC2+BC+BC2+ABC+ABC2+AB2C+AB2C2) > anova(g) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) A 2 34621746 17310873 85.5828 < 2.2e-16 *** B 2 938539 469270 2.3200 0.107992 AB 2 2727451 1363725 6.7421 0.002433 ** AB2 2 570795 285397 1.4110 0.252754 C 2 9549481 4774741 23.6057 4.300e-08 *** AC 2 2985591 1492796 7.3802 0.001467 ** AC2 2 886587 443294 2.1916 0.121580 BC 2 427214 213607 1.0560 0.354901 BC2 2 21134 10567 0.0522 0.949147 ABC 2 4492927 2246464 11.1062 9.119e-05 *** ABC2 2 263016 131508 0.6502 0.525999 AB2C 2 205537 102768 0.5081 0.604500 AB2C2 2 245439 122720 0.6067 0.548815 Residuals 54 10922599 202270 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > > g=lm(y~(A+B+C+D)^4) > anova(g) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) A 2 34621746 17310873 85.5828 < 2.2e-16 *** B 2 938539 469270 2.3200 0.107992 C 2 9549481 4774741 23.6057 4.300e-08 *** D 2 4492927 2246464 11.1062 9.119e-05 *** A:B 4 3298246 824561 4.0765 0.005846 ** A:C 4 3872179 968045 4.7859 0.002231 ** A:D 4 672653 168163 0.8314 0.511145 B:C 2 21134 10567 0.0522 0.949147 B:D 2 205537 102768 0.5081 0.604500 C:D 2 263016 131508 0.6502 0.525999 Residuals 54 10922599 202270 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > > # g=lm(y~A+B+AB+AB2+C+AC+AC2+BC+BC2+D+AD+AD2+BD+BD2+CD+CD2) > g=lm(y~A+B+AB+AB2+C+AC+AC2+BC+BC2+D+AD+BD+CD) # AD2, BD2, CD2 are aliased with BC, AC, AB > anova(g) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) A 2 34621746 17310873 85.5828 < 2.2e-16 *** B 2 938539 469270 2.3200 0.107992 AB 2 2727451 1363725 6.7421 0.002433 ** AB2 2 570795 285397 1.4110 0.252754 C 2 9549481 4774741 23.6057 4.300e-08 *** AC 2 2985591 1492796 7.3802 0.001467 ** AC2 2 886587 443294 2.1916 0.121580 BC 2 427214 213607 1.0560 0.354901 BC2 2 21134 10567 0.0522 0.949147 D 2 4492927 2246464 11.1062 9.119e-05 *** AD 2 245439 122720 0.6067 0.548815 BD 2 205537 102768 0.5081 0.604500 CD 2 263016 131508 0.6502 0.525999 Residuals 54 10922599 202270 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > summary(g) Call: lm(formula = y ~ A + B + AB + AB2 + C + AC + AC2 + BC + BC2 + D + AD + BD + CD) Residuals: Min 1Q Median 3Q Max -1250.7 -124.3 25.0 167.0 824.0 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6223.074 49.972 124.532 < 2e-16 *** A.L 1116.286 86.553 12.897 < 2e-16 *** A.Q -190.244 86.553 -2.198 0.032257 * B.L 178.689 86.553 2.064 0.043788 * B.Q -53.208 86.553 -0.615 0.541304 AB.L 312.960 86.553 3.616 0.000659 *** AB.Q 55.431 86.553 0.640 0.524604 AB2.L -119.056 86.553 -1.376 0.174650 AB2.Q -83.464 86.553 -0.964 0.339191 C.L -589.544 86.553 -6.811 8.34e-09 *** C.Q -78.248 86.553 -0.904 0.369991 AC.L -295.780 86.553 -3.417 0.001209 ** AC.Q -151.959 86.553 -1.756 0.084816 . AC2.L -163.918 86.553 -1.894 0.063609 . AC2.Q 77.250 86.553 0.893 0.376081 BC.L -111.278 86.553 -1.286 0.204050 BC.Q -58.652 86.553 -0.678 0.500896 BC2.L -23.099 86.553 -0.267 0.790583 BC2.Q -15.786 86.553 -0.182 0.855968 D.L -329.852 86.553 -3.811 0.000357 *** D.Q -240.005 86.553 -2.773 0.007610 ** AD.L -18.254 86.553 -0.211 0.833762 AD.Q -93.580 86.553 -1.081 0.284423 BD.L 15.713 86.553 0.182 0.856618 BD.Q -85.823 86.553 -0.992 0.325835 CD.L -7.228 86.553 -0.084 0.933754 CD.Q 98.433 86.553 1.137 0.260453 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 449.7 on 54 degrees of freedom Multiple R-squared: 0.8414, Adjusted R-squared: 0.765 F-statistic: 11.02 on 26 and 54 DF, p-value: 1.375e-13 > par(mfrow=c(1,1)); halfnormal(2*g$coef[-1], label=T, n=7, xlim=c(-.1, 2.6)) > > par(mfrow=c(1,4)) > ylim=c(5000, 7000) > interaction.plot(A, rep(1, length(A)), y, legend=F, ylim=ylim, main="Main Effect A") > interaction.plot(B, rep(1, length(B)), y, legend=F, ylim=ylim, main="Main Effect B") > interaction.plot(C, rep(1, length(C)), y, legend=F, ylim=ylim, main="Main Effect C") > interaction.plot(D, rep(1, length(D)), y, legend=F, ylim=ylim, main="Main Effect D") > > > par(mfrow=c(2,3)) > ylim=c(4500, 7500) > interaction.plot(B, A, y, fixed=T, ylim=ylim) > interaction.plot(C, A, y, fixed=T, ylim=ylim) > interaction.plot(D, A, y, fixed=T, ylim=ylim) > interaction.plot(C, B, y, fixed=T, ylim=ylim) > interaction.plot(D, B, y, fixed=T, ylim=ylim) > interaction.plot(D, C, y, fixed=T, ylim=ylim) > > par(mfrow=c(2,3)) > ylim=c(4500, 7500) > interaction.plot(A, B, y, fixed=T, ylim=ylim) > interaction.plot(A, C, y, fixed=T, ylim=ylim) > interaction.plot(A, D, y, fixed=T, ylim=ylim) > interaction.plot(B, C, y, fixed=T, ylim=ylim) > interaction.plot(B, D, y, fixed=T, ylim=ylim) > interaction.plot(C, D, y, fixed=T, ylim=ylim) > > ## alterative analysis in Section 6.6 > # indicators for D > D01=(D==1) - (D==0) > D02=(D==2) - (D==0) > D12=(D==2) - (D==1) > # standardize > D01=D01/sqrt(2) > D02=D02/sqrt(2) > D12=D12/sqrt(2) > > > g0=lm(y~(A+B+C+D)^2) > g0=lm(y~(A+B+C)^2+(D01+D02+D12)*(A+B+C)) > anova(g0) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) A 2 34621746 17310873 85.5828 < 2.2e-16 *** B 2 938539 469270 2.3200 0.107992 C 2 9549481 4774741 23.6057 4.300e-08 *** D01 1 49747 49747 0.2459 0.621960 D02 1 4443180 4443180 21.9665 1.921e-05 *** A:B 4 3298246 824561 4.0765 0.005846 ** A:C 4 3872179 968045 4.7859 0.002231 ** B:C 4 448348 112087 0.5541 0.696829 A:D01 2 245439 122720 0.6067 0.548815 B:D01 2 205537 102768 0.5081 0.604500 C:D01 2 263016 131508 0.6502 0.525999 Residuals 54 10922599 202270 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > summary(g0) Call: lm(formula = y ~ (A + B + C)^2 + (D01 + D02 + D12) * (A + B + C)) Residuals: Min 1Q Median 3Q Max -1250.7 -124.3 25.0 167.0 824.0 Coefficients: (13 not defined because of singularities) Estimate Std. Error t value Pr(>|t|) (Intercept) 6223.07 49.97 124.532 < 2e-16 *** A.L 1116.29 86.55 12.897 < 2e-16 *** A.Q -190.24 86.55 -2.198 0.032257 * B.L 178.69 86.55 2.064 0.043788 * B.Q -53.21 86.55 -0.615 0.541304 C.L -589.54 86.55 -6.811 8.34e-09 *** C.Q -78.25 86.55 -0.904 0.369991 D01 277.13 99.94 2.773 0.007610 ** D02 -468.42 99.94 -4.687 1.92e-05 *** D12 NA NA NA NA A.L:B.L -343.33 183.61 -1.870 0.066919 . A.Q:B.L 403.12 183.61 2.196 0.032441 * A.L:B.Q 80.25 183.61 0.437 0.663794 A.Q:B.Q 493.67 183.61 2.689 0.009520 ** A.L:C.L 663.83 183.61 3.616 0.000659 *** A.Q:C.L -102.16 183.61 -0.556 0.580237 A.L:C.Q -139.05 183.61 -0.757 0.452165 A.Q:C.Q -221.50 183.61 -1.206 0.232931 B.L:C.L 245.44 183.61 1.337 0.186897 B.Q:C.L -63.12 183.61 -0.344 0.732334 B.L:C.Q -124.90 183.61 -0.680 0.499249 B.Q:C.Q -215.78 183.61 -1.175 0.245066 A.L:D01 -75.89 212.01 -0.358 0.721779 A.Q:D01 220.87 212.01 1.042 0.302158 B.L:D01 -138.44 212.01 -0.653 0.516525 B.Q:D01 162.81 212.01 0.768 0.445867 C.L:D01 135.89 212.01 0.641 0.524266 C.Q:D01 -199.96 212.01 -0.943 0.349813 A.L:D02 NA NA NA NA A.Q:D02 NA NA NA NA B.L:D02 NA NA NA NA B.Q:D02 NA NA NA NA C.L:D02 NA NA NA NA C.Q:D02 NA NA NA NA A.L:D12 NA NA NA NA A.Q:D12 NA NA NA NA B.L:D12 NA NA NA NA B.Q:D12 NA NA NA NA C.L:D12 NA NA NA NA C.Q:D12 NA NA NA NA --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 449.7 on 54 degrees of freedom Multiple R-squared: 0.8414, Adjusted R-squared: 0.765 F-statistic: 11.02 on 26 and 54 DF, p-value: 1.375e-13 > > # stepwise selection using AIC or BIC, the resulting model depends on many choices > # set trace=1 to view the models considered > g=step(g0, direction="both", trace=0, k=2) # AIC > summary(g) Call: lm(formula = y ~ A + B + C + D01 + D02 + A:B + A:C) Residuals: Min 1Q Median 3Q Max -1245.2 -172.2 0.7 226.7 935.8 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6223.074 48.282 128.889 < 2e-16 *** A.L 1116.286 83.628 13.348 < 2e-16 *** A.Q -190.244 83.628 -2.275 0.02627 * B.L 178.689 83.628 2.137 0.03645 * B.Q -53.208 83.628 -0.636 0.52688 C.L -589.544 83.628 -7.050 1.52e-09 *** C.Q -78.248 83.628 -0.936 0.35296 D01 277.133 96.565 2.870 0.00556 ** D02 -468.419 96.565 -4.851 8.22e-06 *** A.L:B.L -290.722 144.847 -2.007 0.04897 * A.Q:B.L 294.288 144.847 2.032 0.04634 * A.L:B.Q -28.579 144.847 -0.197 0.84422 A.Q:B.Q 441.056 144.847 3.045 0.00338 ** A.L:C.L 627.944 144.847 4.335 5.26e-05 *** A.Q:C.L -1.508 144.847 -0.010 0.99173 A.L:C.Q -38.394 144.847 -0.265 0.79181 A.Q:C.Q -185.611 144.847 -1.281 0.20467 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 434.5 on 64 degrees of freedom Multiple R-squared: 0.8245, Adjusted R-squared: 0.7806 F-statistic: 18.79 on 16 and 64 DF, p-value: < 2.2e-16 > > g=step(g0, direction="both", trace=0, k=log(length(y))) # BIC > summary(g) Call: lm(formula = y ~ A + B + C + D01 + D02 + A:B + A:C) Residuals: Min 1Q Median 3Q Max -1245.2 -172.2 0.7 226.7 935.8 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6223.074 48.282 128.889 < 2e-16 *** A.L 1116.286 83.628 13.348 < 2e-16 *** A.Q -190.244 83.628 -2.275 0.02627 * B.L 178.689 83.628 2.137 0.03645 * B.Q -53.208 83.628 -0.636 0.52688 C.L -589.544 83.628 -7.050 1.52e-09 *** C.Q -78.248 83.628 -0.936 0.35296 D01 277.133 96.565 2.870 0.00556 ** D02 -468.419 96.565 -4.851 8.22e-06 *** A.L:B.L -290.722 144.847 -2.007 0.04897 * A.Q:B.L 294.288 144.847 2.032 0.04634 * A.L:B.Q -28.579 144.847 -0.197 0.84422 A.Q:B.Q 441.056 144.847 3.045 0.00338 ** A.L:C.L 627.944 144.847 4.335 5.26e-05 *** A.Q:C.L -1.508 144.847 -0.010 0.99173 A.L:C.Q -38.394 144.847 -0.265 0.79181 A.Q:C.Q -185.611 144.847 -1.281 0.20467 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 434.5 on 64 degrees of freedom Multiple R-squared: 0.8245, Adjusted R-squared: 0.7806 F-statistic: 18.79 on 16 and 64 DF, p-value: < 2.2e-16 > > g1=lm(y~1) > g=step(g1, scope =list(upper = ~(A+B+C)^2+(D01+D02+D12)*(A+B+C), lower = ~1), + direction="both", k=log(nrow(dat)), trace=0) # BIC > summary(g) Call: lm(formula = y ~ A + C + D12 + C:D12 + A:C) Residuals: Min 1Q Median 3Q Max -1443.33 -177.67 50.28 275.33 951.00 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6223.074 51.376 121.127 < 2e-16 *** A.L 1116.286 88.986 12.544 < 2e-16 *** A.Q -190.244 88.986 -2.138 0.036069 * C.L -589.544 88.986 -6.625 6.31e-09 *** C.Q -78.248 88.986 -0.879 0.382279 D12 -372.776 88.986 -4.189 8.14e-05 *** C.L:D12 -486.444 154.129 -3.156 0.002370 ** C.Q:D12 -141.707 154.129 -0.919 0.361085 A.L:C.L 627.944 154.129 4.074 0.000122 *** A.Q:C.L -1.508 154.129 -0.010 0.992224 A.L:C.Q -38.394 154.129 -0.249 0.804022 A.Q:C.Q -185.611 154.129 -1.204 0.232604 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 462.4 on 69 degrees of freedom Multiple R-squared: 0.7858, Adjusted R-squared: 0.7516 F-statistic: 23.01 on 11 and 69 DF, p-value: < 2.2e-16 > > ## avoid effect heredity used in step() > # the final model in WH for flash > g0=lm(y~(A+B+C)^2+(D01+D02+D12)*(A+B+C)) > m=model.matrix(g0)[,-1] > dim(m) [1] 81 39 > dat.m = as.data.frame(m) > g0=lm(y~., dat.m) > g=step(g0, data=dat.m, direction="both", k=log(nrow(dat)), trace=0) # BIC > summary(g) Call: lm(formula = y ~ A.L + A.Q + B.L + C.L + D01 + D02 + `A.L:B.L` + `A.Q:B.L` + `A.Q:B.Q` + `A.L:C.L`, data = dat.m) Residuals: Min 1Q Median 3Q Max -1326.80 -156.56 28.13 201.12 910.53 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6223.07 47.25 131.712 < 2e-16 *** A.L 1116.29 81.84 13.641 < 2e-16 *** A.Q -190.24 81.84 -2.325 0.02299 * B.L 178.69 81.84 2.184 0.03235 * C.L -589.54 81.84 -7.204 5.3e-10 *** D01 277.13 94.50 2.933 0.00454 ** D02 -468.42 94.50 -4.957 4.8e-06 *** `A.L:B.L` -290.72 141.74 -2.051 0.04401 * `A.Q:B.L` 294.29 141.74 2.076 0.04155 * `A.Q:B.Q` 441.06 141.74 3.112 0.00269 ** `A.L:C.L` 627.94 141.74 4.430 3.4e-05 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 425.2 on 70 degrees of freedom Multiple R-squared: 0.8162, Adjusted R-squared: 0.7899 F-statistic: 31.08 on 10 and 70 DF, p-value: < 2.2e-16 > > # the final model in WH for strength > g1=lm(y~1, dat.m) > g=step(g1, scope =list(upper = ~(A.L+A.Q+B.L+B.Q+C.L+C.Q+D01+D02+D12)^2, lower = ~1), + data=dat.m, direction="both", k=log(nrow(dat)), trace=0) # BIC > summary(g) Call: lm(formula = y ~ A.L + C.L + D12 + A.Q + B.L + D01 + A.L:C.L + C.L:D12 + A.Q:B.L, data = dat.m) Residuals: Min 1Q Median 3Q Max -1473.81 -174.98 51.43 237.26 776.27 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6223.07 47.62 130.693 < 2e-16 *** A.L 1116.29 82.47 13.535 < 2e-16 *** C.L -589.54 82.47 -7.148 6.27e-10 *** D12 -468.42 95.23 -4.919 5.44e-06 *** A.Q -190.24 82.47 -2.307 0.02399 * B.L 178.69 82.47 2.167 0.03362 * D01 -191.29 95.23 -2.009 0.04838 * A.L:C.L 627.94 142.85 4.396 3.79e-05 *** C.L:D12 -486.44 142.85 -3.405 0.00109 ** A.Q:B.L 294.29 142.85 2.060 0.04305 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 428.5 on 71 degrees of freedom Multiple R-squared: 0.8106, Adjusted R-squared: 0.7866 F-statistic: 33.77 on 9 and 71 DF, p-value: < 2.2e-16 > > # check the residuals for the final model > par(mfrow=c(2,2)); plot(g,c(1:2,4)); boxcox(g) > > > > ## dispersion model > # > s2 = apply(dat0[,6:8],1,sd) # for strength > #s2 = apply(dat0[,9:11],1,sd) # for flash > y = log(s2) > dat=dat0 > > A=as.ordered(dat[,2]) > B=as.ordered(dat[,3]) > C=as.ordered(dat[,4]) > D=as.ordered(dat[,5]) > > AB=as.ordered((dat[,2]+dat[,3]) %%3) > AB2=as.ordered((dat[,2]+2*dat[,3]) %%3) > AC=as.ordered((dat[,2]+dat[,4]) %%3) > AC2=as.ordered((dat[,2]+2*dat[,4]) %%3) > BC=as.ordered((dat[,3]+dat[,4]) %%3) > BC2=as.ordered((dat[,3]+2*dat[,4]) %%3) > ABC=as.ordered((dat[,2]+dat[,3]+dat[,4]) %%3) > ABC2=as.ordered((dat[,2]+dat[,3]+2*dat[,4]) %%3) > AB2C=as.ordered((dat[,2]+2*dat[,3]+dat[,4]) %%3) > AB2C2=as.ordered((dat[,2]+2*dat[,3]+2*dat[,4]) %%3) > > # AB=CD2, AC=BD2, BC=AD2 > AD=as.ordered((dat[,2]+dat[,5]) %%3) > AD2=as.ordered((dat[,2]+2*dat[,5]) %%3) > BD=as.ordered((dat[,3]+dat[,5]) %%3) > BD2=as.ordered((dat[,3]+2*dat[,5]) %%3) > CD=as.ordered((dat[,4]+dat[,5]) %%3) > CD2=as.ordered((dat[,4]+2*dat[,5]) %%3) # =AB > > > g=lm(y~A+B+AB+AB2+C+AC+AC2+BC+BC2+D+AD+BD+CD) # AD2, BD2, CD2 are aliased with BC, AC, AB > anova(g) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) A 2 7.9483 3.9741 B 2 2.8985 1.4493 AB 2 2.2867 1.1434 AB2 2 5.2834 2.6417 C 2 0.3879 0.1939 AC 2 2.5992 1.2996 AC2 2 0.8457 0.4228 BC 2 0.1391 0.0696 BC2 2 0.1753 0.0876 D 2 1.1774 0.5887 AD 2 1.0314 0.5157 BD 2 0.1844 0.0922 CD 2 0.5955 0.2977 Residuals 0 0.0000 Warning message: In anova.lm(g) : ANOVA F-tests on an essentially perfect fit are unreliable > summary(g) Call: lm(formula = y ~ A + B + AB + AB2 + C + AC + AC2 + BC + BC2 + D + AD + BD + CD) Residuals: ALL 27 residuals are 0: no residual degrees of freedom! Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.44723 NA NA NA A.L -0.85133 NA NA NA A.Q 0.39798 NA NA NA B.L -0.11296 NA NA NA B.Q 0.55615 NA NA NA AB.L -0.50197 NA NA NA AB.Q -0.04595 NA NA NA AB2.L 0.75090 NA NA NA AB2.Q 0.15233 NA NA NA C.L 0.05132 NA NA NA C.Q -0.20115 NA NA NA AC.L 0.38485 NA NA NA AC.Q 0.37509 NA NA NA AC2.L -0.19319 NA NA NA AC2.Q 0.23799 NA NA NA BC.L -0.12277 NA NA NA BC.Q -0.01967 NA NA NA BC2.L 0.09802 NA NA NA BC2.Q -0.09934 NA NA NA D.L 0.34998 NA NA NA D.Q 0.09131 NA NA NA AD.L -0.33133 NA NA NA AD.Q -0.06947 NA NA NA BD.L -0.05502 NA NA NA BD.Q 0.13213 NA NA NA CD.L -0.06662 NA NA NA CD.Q 0.24845 NA NA NA Residual standard error: NaN on 0 degrees of freedom Multiple R-squared: 1, Adjusted R-squared: NaN F-statistic: NaN on 26 and 0 DF, p-value: NA > par(mfrow=c(1,1)); halfnormal(2*g$coef[-1], label=T, n=5, xlim=c(-.1, 2.6)) > > par(mfrow=c(1,4)) > ylim=c(5, 6.2) > interaction.plot(A, rep(1, length(A)), y, legend=F, ylim=ylim, main="Main Effect A") > interaction.plot(B, rep(1, length(B)), y, legend=F, ylim=ylim, main="Main Effect B") > interaction.plot(C, rep(1, length(C)), y, legend=F, ylim=ylim, main="Main Effect C") > interaction.plot(D, rep(1, length(D)), y, legend=F, ylim=ylim, main="Main Effect D") > > > par(mfrow=c(2,3)) > ylim=c(3.8, 6.5) > interaction.plot(B, A, y, fixed=T, ylim=ylim) > interaction.plot(C, A, y, fixed=T, ylim=ylim) > interaction.plot(D, A, y, fixed=T, ylim=ylim) > interaction.plot(C, B, y, fixed=T, ylim=ylim) > interaction.plot(D, B, y, fixed=T, ylim=ylim) > interaction.plot(D, C, y, fixed=T, ylim=ylim) > > par(mfrow=c(2,3)) > ylim=c(3.8, 6.5) > interaction.plot(A, B, y, fixed=T, ylim=ylim) > interaction.plot(A, C, y, fixed=T, ylim=ylim) > interaction.plot(A, D, y, fixed=T, ylim=ylim) > interaction.plot(B, C, y, fixed=T, ylim=ylim) > interaction.plot(B, D, y, fixed=T, ylim=ylim) > interaction.plot(C, D, y, fixed=T, ylim=ylim) >