# Design of Experiments with R
#
# Cano et. al. Six Sigma with R - Statistical Engineering for Process Improvement
# Spinger, UseR! Series, New York, 2012
#
# Chapter 11 - Design of Experiments with R
#
# Walter Ponge-Ferreira
# 09.06.2019
#

# Full 2k Factorial Design
pizzaDesign <- expand.grid(flour = gl(2, 1, labels = c("-","+")),
                           salt = gl(2, 1, labels = c("-","+")),
                           bakPow = gl(2, 1, labels = c("-","+")),
                           score = NA)
pizzaDesign

# Randomize order of experiments
pizzaDesign$ord <- sample(1:8, 8)
pizzaDesign[order(pizzaDesign$ord),]

# Blind trail Results
ss.data.doe1 <- data.frame(repl = rep(1:2, each = 8),
                           rbind(pizzaDesign[,-6], pizzaDesign[,-6]))
ss.data.doe1$score <- c(5.33, 6.99, 4.23, 6.61,
                        2.26, 5.75, 3.26, 6.24,
                        5.70, 7.71, 5.13, 6.76,
                        2.79, 4.57, 2.48, 6.18)

ss.data.doe1

# Average scores
aggregate(score ~ flour + salt + bakPow,
          FUN = mean, data = ss.data.doe1)

boxplot(score ~ flour + salt + bakPow, data = ss.data.doe1)
grid()


# Main effects
(res <- aggregate(score ~ flour,
          FUN = mean, data = ss.data.doe1))
diff(res$score)

boxplot(score ~ flour, data = ss.data.doe1)
grid()

(res <- aggregate(score ~ salt,
          FUN = mean, data = ss.data.doe1))
diff(res$score)

boxplot(score ~ salt, data = ss.data.doe1)
grid()

(res <- aggregate(score ~ bakPow,
                  FUN = mean, data = ss.data.doe1))
diff(res$score)

boxplot(score ~ bakPow, data = ss.data.doe1)
grid()

# Interaction Effects
(res <- aggregate(score ~ flour + salt,
                  FUN = mean, data = ss.data.doe1))
diff(res$score[res$salt == "+"])-diff(res$score[res$salt == "-"])

boxplot(score ~ flour + salt, data = ss.data.doe1)
grid()

(res <- aggregate(score ~ flour + bakPow,
                  FUN = mean, data = ss.data.doe1))
diff(res$score[res$bakPow == "+"])-diff(res$score[res$bakPow == "-"])

boxplot(score ~ flour + bakPow, data = ss.data.doe1)
grid()

(res <- aggregate(score ~ salt + bakPow,
                  FUN = mean, data = ss.data.doe1))
diff(res$score[res$salt == "+"])-diff(res$score[res$salt == "-"])

boxplot(score ~ salt + bakPow, data = ss.data.doe1)
grid()

# Testes de Hipóteses
# H0: mu1 == mu2
# H1: mu1 != mu2

t.test(ss.data.doe1$score[ss.data.doe1$flour == "-"], 
       ss.data.doe1$score[ss.data.doe1$flour == "+"],
       alternative = "two.sided",
       conf.level = 0.99)

t.test(ss.data.doe1$score[ss.data.doe1$salt == "-"], 
       ss.data.doe1$score[ss.data.doe1$salt == "+"],
       alternative = "two.sided",
       conf.level = 0.95)

t.test(ss.data.doe1$score[ss.data.doe1$bakPow == "-"], 
       ss.data.doe1$score[ss.data.doe1$bakPow == "+"],
       alternative = "two.sided",
       conf.level = 0.95)

t.test(ss.data.doe1$score[ss.data.doe1$bakPow == "-"], 
       ss.data.doe1$score[ss.data.doe1$bakPow == "+"],
       alternative = "two.sided",
       conf.level = 0.99)

# ANOVA on Linear Regression
doe.model1 <- lm(score ~ flour + salt + bakPow +
                   flour * salt + flour * bakPow +
                   salt * bakPow + flour * salt * bakPow,
                 data = ss.data.doe1)

summary(doe.model1)

# New Model
doe.model2 <- lm(score ~ flour + bakPow,
                 data = ss.data.doe1)
summary(doe.model2)

coef(doe.model2)

as.numeric(mean(ss.data.doe1$score)) - coef(doe.model2)[2]/2 - coef(doe.model2)[3]/2
coef(doe.model2)[1] + coef(doe.model2)[2]

mean(ss.data.doe1$score) + coef(doe.model2)[2]/2 * (1) + coef(doe.model2)[3]/2 * (-1)

predict(doe.model2)

confint(doe.model2)

plot(c(-1,1), ylim = range(ss.data.doe1$score),
     coef(doe.model1)[1] + c(-1,1) * coef(doe.model1)[2],
     type = "b", pch=16)
abline(h=coef(doe.model1)[1])

plot(c(-1,1), ylim = range(ss.data.doe1$score),
     coef(doe.model1)[1] + c(-1,1) * coef(doe.model1)[3],
     type = "b", pch=16)
abline(h=coef(doe.model1)[1])

plot(c(-1,1), ylim = range(ss.data.doe1$score),
     coef(doe.model1)[1] + c(-1,1) * coef(doe.model1)[4],
     type = "b", pch=16)
abline(h=coef(doe.model1)[1])

plot(c(-1,1), ylim = range(ss.data.doe1$score),
     coef(doe.model2)[1] + c(-1,1) * coef(doe.model2)[2],
     type = "b", pch=16)
abline(h=coef(doe.model2)[1])

plot(c(-1,1), ylim = range(ss.data.doe1$score),
     coef(doe.model2)[1] + c(-1,1) * coef(doe.model2)[3],
     type = "b", pch=16)
abline(h=coef(doe.model2)[1])

prinEf <- data.frame(Factor = rep(c("A_Flour", "C_Backing Powder"), each = 2),
                     Level = rep(c(-1,1), 2),
                     Score = c(aggregate(score ~ flour, FUN = mean, data = ss.data.doe1)[,2],
                               aggregate(score ~ bakPow, FUN = mean, data = ss.data.doe1)[,2]))

library(ggplot2)
p <- ggplot(prinEf,
            aes(x = Level, y = Score)) +
  geom_point() +
  geom_line() +
  scale_x_continuous(breaks = c(-1,1)) +
  facet_grid(. ~ Factor)  +
  geom_abline(intercept = mean(ss.data.doe1$score),
              slope = 0, linetype = "dasked")
plot(p)

intEf <- aggregate(score ~ flour + bakPow,
                   FUN = mean, data = ss.data.doe1)

intEf

p <- ggplot(intEf, aes(x = flour, y = score, color = bakPow )) +
  geom_point() +
  geom_line(aes(group = bakPow)) + 
  geom_abline(intercept = mean(ss.data.doe1$score),
         slope = 0,
         linetype = "dashed",
         color = "black")
plot(p)

p <- ggplot(intEf, aes(x = bakPow, y = score, color = flour )) +
  geom_point() +
  geom_line(aes(group = flour)) + 
  geom_abline(intercept = mean(ss.data.doe1$score),
              slope = 0,
              linetype = "dashed",
              color = "black")
plot(p)

residuals(doe.model2)

par(mfrow=c(2,2)) 
plot(doe.model2)
box("outer")

shapiro.test(residuals(doe.model2))