# R code to accompany the Screening Design video
library(FrF2)
library(rsm)
library(daewr)
library(openxlsx)
library(jtools)
library(interactions)

# Create the fractional factorial design with resolution IV
plan <- FrF2(16, factor.names = c("A", "B", "C", "D", "E", "F"), 
             gen = c("ABC", "BCD"), randomize = FALSE)
summary(plan)
colormap(plan, mod=2)

# All third order interactions are assumed to be negligiable
aliasprint(plan)

# Add center points and save the design to Excel worksheet for data collection
plan2 <- add.center(plan, ncenter=2)

# Read in data with initial experiment results, fit first order interaction model
# Plan points with experiment results in file "data/Tab17-2.xlsx"
dat.partial <- read.xlsx("data/Tab17-2.xlsx", sheet="partial")
model.init <- rsm(Loge_Etch_Uniformity~FO(A, B, C, D, E, F)+TWI(A, B, C, D, E, F), data=dat.partial)
# or
# model.init <- lm(Loge_Etch_Uniformity~(A+B+C+D+E+F)^2, data=dat.partial)
summary(model.init)

# Semi-fold to expand the model, collect data, read in the full data worksheet
dat <- read.xlsx("data/Tab17-2.xlsx", sheet="full")
model.full <- rsm(Loge_Etch_Uniformity~FO(A, B, C, D, E, F)+TWI(A, B, C, D, E, F), data=dat)
# The same as
# model.full <- lm(Loge_Etch_Uniformity~(A+B+C+D+E+F)^2+Block, data=dat)
summary(model.full)

## Refind the full model to only those are significant, or the book models. 
# Use the rsm() so we can see more carnonical summary statistics and also contour plots later
model.book <- rsm(Loge_Etch_Uniformity~FO(A, B, C, D, E, F) + Block + B*C+A*F+E*F, data=dat)
summary(model.book)

# The model without the Block, for optimization blocks makes no sense, contours no sense neither
model.book.2 <- rsm(Loge_Etch_Uniformity~FO(A, B, C, D, E, F) + B*C+A*F+E*F, data=dat)
summary(model.book.2)

#############################################################################
## Main effect plots of the final model at the optimal points
effect_plot(model.book.2, pred = E, interval = TRUE, plot.points = TRUE,
            at=summary(model.2)$canonical$xs,
            x.label = "E: 灼蚀气体流速", y.label = "氧化硅灼蚀一致性")
effect_plot(model.book.2, pred = F, interval = TRUE, plot.points = TRUE,
            at=summary(model.2)$canonical$xs,
            x.label = "F: 氧化硅灼蚀厚度 (50 A vs. 200 A)", y.label = " 氧化硅灼蚀一致性")

## Interaction effect plots of final model at the optimal points
interact_plot(model.book.2, pred = B, modx = C, modx.values=c(-1, 1), 
              interval = TRUE, plot.points = TRUE,
              at=summary(model.book.2)$canonical$xs,
              x.label = "B: 灼蚀前氮气和氮水混合气总气流量", 
              modx.labels = c("-1: 低流量", "+1: 高流量"),
              y.label="氧化硅灼蚀一致性",
              legend.main="C: 灼蚀前水气雾流量")
interact_plot(model.book.2, pred = A, modx = F, modx.values=c(-1, 1), 
              interval = TRUE, plot.points = TRUE,
              at=summary(model.book.2)$canonical$xs,
              x.label = "A: 托盘转速", 
              modx.labels = c("-1: 50 Ang.", "+1: 200 Ang."),
              y.label="氧化硅灼蚀一致性",
              legend.main="F: 氧化硅灼蚀厚度 \n(50 A vs. 200 A)")
interact_plot(model.book.2, pred = E, modx = F, modx.values=c(-1, 1), 
              interval = TRUE, plot.points = TRUE,
              at=summary(model.book.2)$canonical$xs,
              x.label = "E: 灼蚀气体流速", 
              modx.labels = c("-1: 50 Ang.", "+1: 200 Ang."),
              y.label="氧化硅灼蚀一致性",
              legend.main="F: 氧化硅灼蚀厚度 \n(50 A vs. 200 A)")

#############################################################################
# Contour plots at the optimal condition
# Figure 17.3, A: 托盘转速 vs. F: 氧化硅灼蚀厚度 (50 A vs. 200 A)
par(mfrow=c(1,2))
contour(model.book.2, ~ A+F, image = TRUE, at = summary(model.2)$canonical$xs)
persp(model.book.2, ~ A+F, 
      contour=list(Loge_Etch_Uniformity="bottom"), at = summary(model.2)$canonical$xs, 
      theta=300,
      phi=25)

# Figure 17.4, C: 灼蚀前水气雾流量 vs. B: 灼蚀前氮气和氮水混合气总气流量
par(mfrow=c(1,2))
contour(model.book.2, ~ C+B, image = TRUE, at = summary(model.2)$canonical$xs)
persp(model.book.2, ~ C+B, 
      contour=list(Loge_Etch_Uniformity="bottom"), at = summary(model.2)$canonical$xs, 
      theta=300,
      phi=25)

# Figure 17.5, F:氧化硅灼蚀厚度 (50 A vs. 200 A) vs. E:灼蚀气体流速,
par(mfrow=c(1,2))
contour(model.book.2, ~ F+E, image = TRUE, at = summary(model.2)$canonical$xs, nlevels=20)
persp(model.book.2, ~ F+E, 
      contour=list(Loge_Etch_Uniformity="bottom"), at = summary(model.2)$canonical$xs, 
      theta=300,
      phi=25)