pdf("C:\\Rmisc\\graphs\\orange_emul.pdf") orange.dat <- read.fwf("C:\\data\\orange_emul.txt",width=rep(8,12), col.names=c("block","x1","x2","x3","block1","block2","block3","y1","y2","y3","y4","y5")) attach(orange.dat) orange.dat x11 <- x1^2; x22 <- x2^2; x33 <- x3^2; x12 <- x1*x2; x13 <- x1*x3; x23 <- x2*x3; mean(y2) orange.rsm1 <- lm(y2 ~ x1+x2+x3+x11+x22+x33+x12+x13+x23+block1+block2) summary(orange.rsm1) orange.rsm2 <- lm(y2 ~ x1+x2+x3+x11+x22+x33+x12+x13+x23) summary(orange.rsm2) anova(orange.rsm1,orange.rsm2) x1.lo <- 13; x1.hi <- 20; x1.step <- 0.05; x1.rows <- 1+(x1.hi-x1.lo)/x1.step; x2.lo <- 0.3; x2.hi <- 0.5; x2.step <- 0.01; x2.cols <- 1+(x2.hi-x2.lo)/x2.step; x3.lo <- 10; x3.hi <- 14; x3.step <- 0.04; x3.cols <- 1+(x3.hi-x3.lo)/x3.step; ########### Plot Response surface versus (X1,X2) at various X3 Values yhatx1x2.1 <- matrix(rep(0,x1.rows*x2.cols),ncol=x2.cols) yhatx1x2.2 <- matrix(rep(0,x1.rows*x2.cols),ncol=x2.cols) yhatx1x2.3 <- matrix(rep(0,x1.rows*x2.cols),ncol=x2.cols) x1.elem <- 0 x2.elem <- 0 x3.1 <- 10; x3.2 <- 12; x3.3 <- 14; for (x1 in seq(x1.lo,x1.hi,x1.step)) { x1.elem <- x1.elem+1 x2.elem <- 0 for (x2 in seq(x2.lo,x2.hi,x2.step)) { x2.elem <- x2.elem+1 yhatx1x2.1[x1.elem,x2.elem] <- (80.4006+49.5912*x1-2106.4948*x2+119.0082*x3.1+ 0.1287*(x1^2)+480.4509*(x2^2)-3.3939*(x3.1^2)+ 37.50*x1*x2-3.7679*x1*x3.1+150.6250*x2*x3.1) yhatx1x2.2[x1.elem,x2.elem] <- (80.4006+49.5912*x1-2106.4948*x2+119.0082*x3.2+ 0.1287*(x1^2)+480.4509*(x2^2)-3.3939*(x3.2^2)+ 37.50*x1*x2-3.7679*x1*x3.2+150.6250*x2*x3.2) yhatx1x2.3[x1.elem,x2.elem] <- (80.4006+49.5912*x1-2106.4948*x2+119.0082*x3.3+ 0.1287*(x1^2)+480.4509*(x2^2)-3.3939*(x3.3^2)+ 37.50*x1*x2-3.7679*x1*x3.3+150.6250*x2*x3.3) }} x1 <- seq(x1.lo,x1.hi,x1.step) x2 <- seq(x2.lo,x2.hi,x2.step) persp(x1,x2,yhatx1x2.1,xlim=c(13,20),ylim=c(0.3,0.5),main="y-hat vs x1,x2 @ x3=10") persp(x1,x2,yhatx1x2.2,xlim=c(13,20),ylim=c(0.3,0.5),main="y-hat vs x1,x2 @ x3=12") persp(x1,x2,yhatx1x2.3,xlim=c(13,20),ylim=c(0.3,0.5),main="y-hat vs x1,x2 @ x3=14") contour(x1,x2,yhatx1x2.1,xlim=c(13,20),ylim=c(0.3,0.5),main="y-hat vs x1,x2 @ x3=10") contour(x1,x2,yhatx1x2.2,xlim=c(13,20),ylim=c(0.3,0.5),main="y-hat vs x1,x2 @ x3=12") contour(x1,x2,yhatx1x2.3,xlim=c(13,20),ylim=c(0.3,0.5),main="y-hat vs x1,x2 @ x3=14") ########### Plot Response surface versus (X1,X3) at various X2 Values yhatx1x3.1 <- matrix(rep(0,x1.rows*x3.cols),ncol=x3.cols) yhatx1x3.2 <- matrix(rep(0,x1.rows*x3.cols),ncol=x3.cols) yhatx1x3.3 <- matrix(rep(0,x1.rows*x3.cols),ncol=x3.cols) x1.elem <- 0 x3.elem <- 0 x2.1 <- 0.3; x2.2 <- 0.4; x2.3 <- 0.5; for (x1 in seq(x1.lo,x1.hi,x1.step)) { x1.elem <- x1.elem+1 x3.elem <- 0 for (x3 in seq(x3.lo,x3.hi,x3.step)) { x3.elem <- x3.elem+1 yhatx1x3.1[x1.elem,x3.elem] <- (80.4006+49.5912*x1-2106.4948*x2.1+119.0082*x3+ 0.1287*(x1^2)+480.4509*(x2.1^2)-3.3939*(x3^2)+ 37.50*x1*x2.1-3.7679*x1*x3+150.6250*x2.1*x3) yhatx1x3.2[x1.elem,x3.elem] <- (80.4006+49.5912*x1-2106.4948*x2.2+119.0082*x3+ 0.1287*(x1^2)+480.4509*(x2.2^2)-3.3939*(x3^2)+ 37.50*x1*x2.2-3.7679*x1*x3+150.6250*x2.2*x3) yhatx1x3.3[x1.elem,x3.elem] <- (80.4006+49.5912*x1-2106.4948*x2.3+119.0082*x3+ 0.1287*(x1^2)+480.4509*(x2.3^2)-3.3939*(x3^2)+ 37.50*x1*x2.3-3.7679*x1*x3+150.6250*x2.3*x3) }} x1 <- seq(x1.lo,x1.hi,x1.step) x3 <- seq(x3.lo,x3.hi,x3.step) persp(x1,x3,yhatx1x3.1,xlim=c(13,20),ylim=c(10,14),main="y-hat vs x1,x3 @ x2=0.3") persp(x1,x3,yhatx1x3.2,xlim=c(13,20),ylim=c(10,14),main="y-hat vs x1,x3 @ x2=0.4") persp(x1,x3,yhatx1x3.3,xlim=c(13,20),ylim=c(10,14),main="y-hat vs x1,x3 @ x2=0.5") contour(x1,x3,yhatx1x3.1,xlim=c(13,20),ylim=c(10,14),main="y-hat vs x1,x3 @ x2=0.3") contour(x1,x3,yhatx1x3.2,xlim=c(13,20),ylim=c(10,14),main="y-hat vs x1,x3 @ x2=0.4") contour(x1,x3,yhatx1x3.3,xlim=c(13,20),ylim=c(10,14),main="y-hat vs x1,x3 @ x2=0.5") ########### Plot Response surface versus (X2,X3) at various X1 Values x2.rows=x2.cols yhatx2x3.1 <- matrix(rep(0,x2.rows*x3.cols),ncol=x3.cols) yhatx2x3.2 <- matrix(rep(0,x2.rows*x3.cols),ncol=x3.cols) yhatx2x3.3 <- matrix(rep(0,x2.rows*x3.cols),ncol=x3.cols) x2.elem <- 0 x3.elem <- 0 x1.1 <- 13; x1.2 <- 16.5; x1.3 <- 20; for (x2 in seq(x2.lo,x2.hi,x2.step)) { x2.elem <- x2.elem+1 x3.elem <- 0 for (x3 in seq(x3.lo,x3.hi,x3.step)) { x3.elem <- x3.elem+1 yhatx2x3.1[x2.elem,x3.elem] <- (80.4006+49.5912*x1.1-2106.4948*x2+119.0082*x3+ 0.1287*(x1.1^2)+480.4509*(x2^2)-3.3939*(x3^2)+ 37.50*x1.1*x2-3.7679*x1.1*x3+150.6250*x2*x3) yhatx2x3.2[x2.elem,x3.elem] <- (80.4006+49.5912*x1.2-2106.4948*x2+119.0082*x3+ 0.1287*(x1.2^2)+480.4509*(x2^2)-3.3939*(x3^2)+ 37.50*x1.2*x2-3.7679*x1.2*x3+150.6250*x2*x3) yhatx2x3.3[x2.elem,x3.elem] <- (80.4006+49.5912*x1.3-2106.4948*x2+119.0082*x3+ 0.1287*(x1.3^2)+480.4509*(x2^2)-3.3939*(x3^2)+ 37.50*x1.3*x2-3.7679*x1.3*x3+150.6250*x2*x3) }} x2 <- seq(x2.lo,x2.hi,x2.step) x3 <- seq(x3.lo,x3.hi,x3.step) persp(x2,x3,yhatx2x3.1,xlim=c(0.3,0.5),ylim=c(10,14),main="y-hat vs x2,x3 @ x1=13") persp(x2,x3,yhatx2x3.2,xlim=c(0.3,0.5),ylim=c(10,14),main="y-hat vs x2,x3 @ x1=16.5") persp(x2,x3,yhatx2x3.3,xlim=c(0.3,0.5),ylim=c(10,14),main="y-hat vs x2,x3 @ x1=20") contour(x2,x3,yhatx2x3.1,xlim=c(0.3,0.5),ylim=c(10,14),main="y-hat vs x2,x3 @ x1=13") contour(x2,x3,yhatx2x3.2,xlim=c(0.3,0.5),ylim=c(10,14),main="y-hat vs x2,x3 @ x1=16.5") contour(x2,x3,yhatx2x3.3,xlim=c(0.3,0.5),ylim=c(10,14),main="y-hat vs x2,x3 @ x1=20") dev.off()