# Problem 14.1 (p. 515)
impure=read.table("impurity.txt",header=T)
attach(impure)
# fit response surface using all two-way interactions and all quadratic terms
# include a block effect here -- usually don't block in computer experiments
impure.lm=lm(Impurity~X1*X2*X3+I(X1^2)+I(X2^2)+I(X3^2)+Block)
summary(impure.lm)

      Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
      (Intercept)  18.9614     0.4876  38.891 2.10e-10 ***
      X1            6.1787     0.3234  19.105 5.84e-08 ***
      X2           -0.1047     0.3234  -0.324   0.7545    
      X3            4.0629     0.3234  12.563 1.51e-06 ***
      I(X1^2)       0.7234     0.3160   2.289   0.0513 .  
      I(X2^2)       0.3992     0.3160   1.263   0.2420    
      I(X3^2)       4.3373     0.3160  13.727 7.65e-07 ***
      Block         0.3441     0.2702   1.274   0.2385    
      X1:X2         0.2412     0.4224   0.571   0.5836    
      X1:X3        -5.0863     0.4224 -12.042 2.09e-06 ***
      X2:X3         0.2462     0.4224   0.583   0.5759    
      X1:X2:X3      0.2737     0.4224   0.648   0.5351    

# lots of insignificant terms, but first we should check assumptions
# check independence assumption by plotting residuals vs. run order
plot(RunOrder,resid(impure.lm))  # looks fine
# check quadratic/linear fit, and check equal variance assumption
pairs(cbind(resid(impure.lm),impure[,3:5]))
plot(fitted(impure.lm),resid(impure.lm))
# hard to tell in a central composite design; this works better for 3^k designs
# check normality assumption
qqnorm(resid(impure.lm))    # look OK

# since assumptions look OK, now do model selection, removing terms
#  one-at-a-time and refitting the model after each drop
#  keep intercept and keep main effects if their interactions or
#  quadratic terms are kept
# so the first term we drop is the 3-way interaction
impure.lm2=lm(Impurity~(X1+X2+X3)^2+I(X1^2)+I(X2^2)+I(X3^2)+Block)
summary(impure.lm2)
# now drop X1:X2
impure.lm3=lm(Impurity~(X1+X2+X3)^2-X1:X2+I(X1^2)+I(X2^2)+I(X3^2)+Block)
summary(impure.lm3)
# now drop X2:X3
impure.lm4=lm(Impurity~X1+X2+X3+X1:X3+I(X1^2)+I(X2^2)+I(X3^2)+Block)
summary(impure.lm4)
# now drop X2^2
impure.lm5=lm(Impurity~X1+X2+X3+X1:X3+I(X1^2)+I(X3^2)+Block)
summary(impure.lm5)
# finally we can drop the main effect for X2
impure.lm6=lm(Impurity~X1+X3+X1:X3+I(X1^2)+I(X3^2)+Block)
summary(impure.lm6)
# and now we can drop the blocking effect
impure.lm7=lm(Impurity~X1+X3+X1:X3+I(X1^2)+I(X3^2))
summary(impure.lm7)
# double check residuals
qqnorm(resid(impure.lm7))    # still look OK


# plotting the surface
# first make a grid for X1 and X3 over which we'll get fitted values
impure.grid=expand.grid(X1=seq(-1.75,1.75,by=.25),X3=seq(-1.75,1.75,by=.25))
# get fitted values over this grid
impure.predict=predict(impure.lm7,impure.grid)
# 3-D plots require the X and Y arguments to be sequences of the grid values
#  (but just the unique sorted values), and the Z to be a matrix
contour(seq(-1.75,1.75,by=.25),seq(-1.75,1.75,by=.25),matrix(impure.predict,nrow=15,ncol=15))
persp(seq(-1.75,1.75,by=.25),seq(-1.75,1.75,by=.25),matrix(impure.predict,nrow=15,ncol=15))

# finding an optimum
# first define a function in R for the fitted surface
#  function must be of a single vector
coefficients(impure.lm7)
impure.function=function(x) {
  19.2878682 + 6.1787348*x[1] + 4.0628811*x[2] + 0.6585319*x[1]^2 + 4.2724775*x[2]^2 - 5.0862500*x[1]*x[2]
}
# now run "non-linear minimization
#  (if we wanted a maximum, use the negative of the function in nlm)
#  arguments are the function and the initial values for the search
nlm(impure.function,c(0,0))
# that fails, but we can see that we want small X1 and small X3,
#  so over the search region, we'll call the minimum (-1.75, -1.75)

# 95% confidence interval for fitted value at minimum
predict(impure.lm7,data.frame(X1=-1.75,X3=-1.75),interval="confidence")
# 95% confidence interval for predicted value at minimum
predict(impure.lm7,data.frame(X1=-1.75,X3=-1.75),interval="predict")