# Analysis of a Central Composite Design for 2 Factors

# Author: Rosane Rech, January 2021.
# This code is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
# https://creativecommons.org/licenses/by-nc-sa/4.0/ 

# Data source:
# Design and analysis of experiments / Douglas C. Montgomery. — Eighth edition.
# ISBN 978-1-118-14692-7

# Data file: DoEOpt06
# Time (x1): Time (min)
# Temp (x2): Temperature (ºC)
# Y: Yield (%)


# building the data set 
# (run this code to build the DoEOpt05 data set before following the analysis)

DoEOpt06 <- data.frame(x1 = c(-1, -1, 1, 1, 0, 0, 0, 0, 0, 1.414, -1.414, 0, 0),
                          x2 = c(-1, 1, -1, 1, 0, 0, 0, 0, 0, 0, 0, 1.414, -1.414),
                          Time = c(80, 80, 90, 90, 85, 85, 85, 85, 85, 92.07, 77.93, 85, 85),
                          Temp = c(170, 180, 170, 180, 175, 175, 175, 175, 175, 175, 175, 182.07, 167.93),
                          Y = c(76.5, 77, 78, 79.5, 79.9, 80.3, 80, 79.7, 79.8, 78.4, 75.6, 78.5, 77)
                            )

# loading Response Surface Methodology package

library(rsm)

# checking file structure

str(DoEOpt06)

# setting the realtionship between the coded and the natural variables

DoEOpt06 <- as.coded.data(DoEOpt06, 
              x1 ~ (Time-85)/5,
              x2 ~ (Temp-175)/5)

###
##
# regression model for the Yield

model_Y <- rsm(Y ~ SO(x1,x2), data = DoEOpt06)

model_Y <- rsm(Y ~ FO(x1,x2) + TWI(x1,x2) + PQ(x1,x2), data = DoEOpt06)
summary(model_Y)

# contour and perspective plots

contour(model_Y, x1~x2, image = TRUE, 
        xlabs=c("Time (min)", "Temperature (ºC)"))

persp(model_Y, x1~x2, col = terrain.colors(50), contours = "colors",
      zlab = "Yield (%)", 
      xlabs=c("Time (min)", "Temperature (ºC)"))

# predictig the Yield at the stationary point

max <- data.frame(x1 = 0.361, x2 = 0.257)
predict(model_Y, max)