## Chemical process Experiment y=c(28, 25, 27, 36, 32, 32 , 18, 19, 23, 31, 30, 29) A=rep(rep(c(-1,1), c(3,3)), 2) B=rep(rep(c(-1,1),c(6,6)), 1) g=lm(y~A*B) summary(g) anova(g) model.matrix(g) g=lm(y~A+B) summary(g) anova(g) par(mfrow=c(3,2)) interaction.plot(A, rep(1, length(A)), y, legend=F, main="Main Effect A") interaction.plot(B, rep(1, length(B)), y, legend=F, main="Main Effect B") interaction.plot(A, B, y, fixed=T, main="Interaction AB") interaction.plot(B, A, y, fixed=T, main="Interaction BA") plot(g, 1:2) ## if there are missing observations so that the design is unbalanced y[12] = NA # set last obs missing g1=lm(y~A*B) summary(g1) anova(g1) g2=lm(y~B*A) anova(g2) ## Filtration Rate Experiment y=c(45, 71, 48, 65, 68, 60, 80, 65, 43, 100, 45, 104, 75, 86, 70, 96) A=rep(c(-1,1), 8) B=rep(rep(c(-1,1), c(2,2)), 4) C=rep(rep(c(-1,1), c(4,4)), 2) D=rep(c(-1,1), c(8,8)) g=lm(y~(A+B+C+D)^4) summary(g) source("halfnormal.R") # download from Web and save it in working directory source("http://www.stat.ucla.edu/~hqxu/stat201A/R/halfnormal.R") # if online par(mfrow=c(1,2)) normalplot(2*coef(g)[-1], l=T, n=5, ylim=c(-20, 25)) halfnormalplot(2*coef(g)[-1], l=T, n=5) g2=lm(y~(A+C+D)^3) summary(g2) g3=lm(y~A+C+D+A:C +A:D) summary(g3) # residual plots for the final model par(mfrow=c(1,2)) plot(g3, 1:2) ## main effect and interaction plots par(mfrow=c(3,3)) interaction.plot(A, rep(1, length(y)), y, ylim=c(40, 100),legend=F, main="Main Effect A") interaction.plot(C, rep(1, length(y)), y, ylim=c(40, 100),legend=F, main="Main Effect C") interaction.plot(D, rep(1, length(y)), y, ylim=c(40, 100),legend=F, main="Main Effect D") interaction.plot(A, C, y, ylim=c(40, 100), fixed=T, main="Interaction AC") interaction.plot(A, D, y, ylim=c(40, 100), fixed=T, main="Interaction AD") interaction.plot(C, D, y, ylim=c(40, 100), fixed=T, main="Interaction CD") CD=as.factor(paste(C, D, sep="")) interaction.plot(A, CD, y, ylim=range(y), fixed=T, main="Interaction ACD") AC=as.factor(paste(A, C, sep="")) interaction.plot(D, AC, y, ylim=range(y), fixed=T, main="Interaction ACD") AD=as.factor(paste(A, D, sep="")) interaction.plot(C, AD, y, ylim=range(y), fixed=T, main="Interaction ACD") ## Chemical process Experiment in 3 blocks y=c(28, 25, 27, 36, 32, 32 , 18, 19, 23, 31, 30, 29) A=rep(rep(c(-1,1), c(3,3)), 2) B=rep(rep(c(-1,1),c(6,6)), 1) block=as.factor(rep(c(1,2,3), 4)) g=lm(y~block+A*B) anova(g) ## Filtration Rate Experiment in 2 blocks # y=c(45, 71, 48, 65, 68, 60, 80, 65, 43, 100, 45, 104, 75, 86, 70, 96) y=c(25,71,48,45,68,40,60,65,43,80,25,104,55,86,70,76) # block1 -20 block=as.factor(c(1,2,2,1,2,1,1,2,2,1,1,2,1,2,2,1)) A=rep(c(-1,1), 8) B=rep(rep(c(-1,1), c(2,2)), 4) C=rep(rep(c(-1,1), c(4,4)), 2) D=rep(c(-1,1), c(8,8)) options(contrasts=c("contr.sum", "contr.poly")) # zero-sum constraint for block g=lm(y~(A+B+C+D)^3+block) summary(g) 2*coef(g) ## A half-fraction of the Filtration Rate Experiment, I=ABCD y=c(45, 71, 48, 65, 68, 60, 80, 65, 43, 100, 45, 104, 75, 86, 70, 96) A=rep(c(-1,1), 8) B=rep(rep(c(-1,1), c(2,2)), 4) C=rep(rep(c(-1,1), c(4,4)), 2) D=rep(c(-1,1), c(8,8)) dat=data.frame(cbind(A,B,C,D,y)) subset= D==A*B*C # half-fraction, I=ABCD dat[subset,] g=lm(y~A+B+C+D+A:B+A:C+A:D, dat, subset=subset) summary(g) 2*coef(g)