# first, we need to enter the data used in SPS06, which we do
# using coded values for the factors

temp = c(1, 1, 1, 1, -1, -1, -1, -1)

ppd = c(-1, 1, 1, -1, -1, 1, 1, -1)

ph = c(1, 1, -1, -1, 1, 1, -1, -1)

rate = c(6.69, 11.71, 14.79, 8.05, 6.33, 11.11, 14.08, 7.59)

# second, we build a model of the response surface; note that the
# use of * means that the model will include first-order effects
# for each facator and all possible interactions between the
# factors; because the number of parameters in the model (8) is
# the same as the number of data points, a summary yields not
# useful statistical information

cer.lm = lm(rate ~ temp * ppd * ph)

summary(cer.lm)

# if our data is in a data frame, then we need to tell R where to
# find the data, which we can do using data = data.frame or by
# using data.frame$column name
cer = data.frame(temp, ppd, ph, rate)

cer.lm2 = lm(rate ~ temp * ppd * ph, data = cer)

summary(cer.lm2)

cer.lm3 = lm(cer$rate ~ cer$temp * cer$ppd * cer$ph)

summary(cer.lm3)

# third, we evaluate the model to identify important
# coefficients; if the responses are just random noise
# independent of the factors, then we expect that the
# coefficients should be small in magnitude and distributed
# randomly (that is, a normal distribution) around a mean given
# by the y-intercept; thus, they should fall along a qqline on a
# qqplot

qqnorm(coef(cer.lm)[-1])

qqline(coef(cer.lm)[-1])

# fourth, we refine the model by keeping only those parameters
# that seem likely to have a real affect on the response; in this
# case that is ppd, ph, and their interaction

cer.new = lm(rate ~ ppd * ph)

summary(cer.new)

# one limitation to our approach is that by testing each factor
# at only two levels, we are limited to a model that is a
# straight-line; if we wish to evalute if this is reasonable,
# then we can make replicate measurements at the center of our
# experimental design and use a t-test to see if there is a
# significant difference between the y-intercept and the mean of
# the replciates finding, in this case, that there is no evidence
# of a significant difference

cer.rep = c(10.06, 9.81, 10.04, 10.16, 10.24)

mean(cer.rep)

sd(cer.rep)

t.test(cer.rep, mu = coef(cer.new)[1], alternative = "two.sided")

# fifth, once we have a model we can use it to predict the
# response for combinations of factor levels; we do this using
# the predict function

help(predict.lm)

cer.pred = data.frame(temp = -0.600, ppd = -0.291, ph = 0.75)

predict(cer.new, cer.pred)

# sixth, visualizing a model is always a good idea; a plot of the
# response versus two factors requires three-dimensions with the
# response on the z-axis and the factors on the x-axis and the
# y-axis; we can do this using the plot3D package; first we
# create vectors to hold values for ppd and for ph that we can
# use to calculate the response

ppd.x = seq(-1, 1, 0.1)

ph.y = seq(-1, 1, 0.1)

# next, we write a function that will take a value for x anf for
# y and return the response

model = function(x, y) {coef(cer.new)[1] + coef(cer.new)[2]*x + 
    coef(cer.new)[3]*y + coef(cer.new)[4]*x*y}

# finally, we use the outer function to calculate the response
# for every possible combination of x and of y
help(outer)

rate.z = outer(ppd.x, ph.y, model)

View(rate.z)

# now we can plot the response surface using perspective plots or contour plots

library(plot3D)

persp3D(x = ppd.x, y = ph.y, z = rate.z, phi = 0, theta = 30)

persp3D(x = ppd.x, y = ph.y, z = rate.z, phi = 30, theta = 30)

persp3D(x = ppd.x, y = ph.y, z = rate.z, phi = 0, theta = 60)

contour2D(x = ppd.x, y = ph.y, z = rate.z)

persp3D(x = ppd.x, y = ph.y, z = rate.z, phi = 0, theta = 30, contour = TRUE, ticktype = "detailed")