Do the following problem in R
In a process development study on yield, four factors were studied, each at two levels: time(A), concentration(B), pressure(C), and temperature(D).
# Read in the data table for 14.1
dat16<-read.csv("https://raw.githubusercontent.com/forestwhite/RStatistics/main/16table.csv")
dat16[dat16 == "+"] <- 1
dat16[dat16 == "-"] <- -1
dat16$A<- as.numeric(dat16$A)
dat16$B<- as.numeric(dat16$B)
dat16$C<- as.numeric(dat16$C)
dat16$D<- as.numeric(dat16$D)
dat16
## RunNumber ActualRunOrder A B C D Yield
## 1 1 5 -1 -1 -1 -1 12
## 2 2 9 1 -1 -1 -1 18
## 3 3 8 -1 1 -1 -1 13
## 4 4 13 1 1 -1 -1 16
## 5 5 3 -1 -1 1 -1 17
## 6 6 7 1 -1 1 -1 15
## 7 7 14 -1 1 1 -1 20
## 8 8 1 1 1 1 -1 15
## 9 9 6 -1 -1 -1 1 10
## 10 10 11 1 -1 -1 1 25
## 11 11 2 -1 1 -1 1 13
## 12 12 15 1 1 -1 1 24
## 13 13 4 -1 -1 1 1 19
## 14 14 16 1 -1 1 1 21
## 15 15 10 -1 1 1 1 17
## 16 16 12 1 1 1 1 23
(a) Display the halfnormal plot for this data and determine which factors appear to be significant.
Significant effects that differ from the norm are time(A), temperature(D), the interaction of time(A) and pressure(C), and the interaction of time(A) and temperature(D).
library(DoE.base)
## Loading required package: grid
## Loading required package: conf.design
## Registered S3 method overwritten by 'DoE.base':
## method from
## factorize.factor conf.design
##
## Attaching package: 'DoE.base'
## The following objects are masked from 'package:stats':
##
## aov, lm
## The following object is masked from 'package:graphics':
##
## plot.design
## The following object is masked from 'package:base':
##
## lengths
mod16<-lm(Yield~A*B*C*D,data=dat16)
coef(mod16)
## (Intercept) A B C D
## 1.737500e+01 2.250000e+00 2.500000e-01 1.000000e+00 1.625000e+00
## A:B A:C B:C A:D B:D
## -3.750000e-01 -2.125000e+00 1.250000e-01 2.000000e+00 -8.543513e-17
## C:D A:B:C A:B:D A:C:D B:C:D
## 1.110223e-16 5.000000e-01 3.750000e-01 -1.250000e-01 -3.750000e-01
## A:B:C:D
## 5.000000e-01
halfnormal(mod16, ME.partial=FALSE)
##
## Significant effects (alpha=0.05, Lenth method):
## [1] A A:C A:D D
(b) Pull terms that do not appear to be significant into error and test for the significance of the other effects at the 0.05 level of significance.
Pulling out concentration(B) entirely, pressure(C) on its own,and interaction effects for pressure(C) and temperature(D), we can reject H0 for time(A), temperature(C), temperature(D), the interaction of time(A) and pressure(C), and the interaction of time(A) and temperature(D). Analysis of the residuals indicates that the model is viable (constant variance, normality). The interaction plots indicate that high (1) time has a positive correlation with Yield, but pressure and temperature interfere with the time effect.
# designate fixed factors as fixed
library(GAD)
## Loading required package: matrixStats
## Loading required package: R.methodsS3
## R.methodsS3 v1.8.1 (2020-08-26 16:20:06 UTC) successfully loaded. See ?R.methodsS3 for help.
Time <- as.fixed(dat16$A)
Time
## [1] -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1
## Levels: -1 1
Pressure <- as.fixed(dat16$C)
Pressure
## [1] -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 1 1
## Levels: -1 1
Temperature <- as.fixed(dat16$D)
Temperature
## [1] -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 1 1 1
## Levels: -1 1
# Apply analysis of variance model
model16 <- lm(dat16$Yield ~ Time+Temperature+Time*Pressure+Time*Temperature)
gad(model16)
## Analysis of Variance Table
##
## Response: dat16$Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Time 1 81.00 81.000 49.8462 3.456e-05 ***
## Temperature 1 42.25 42.250 26.0000 0.0004647 ***
## Pressure 1 16.00 16.000 9.8462 0.0105485 *
## Time:Pressure 1 72.25 72.250 44.4615 5.583e-05 ***
## Time:Temperature 1 64.00 64.000 39.3846 9.193e-05 ***
## Residual 10 16.25 1.625
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(model16,1)
plot(model16,2)
# Produce an interaction plot for Time*Pressure
interaction.plot(x.factor = Time,
trace.factor = Pressure,
response = dat16$Yield,
fun = mean,
type="b", ### Show plot points as symbols
col=c("steelblue","firebrick2"), ### Colors for levels of trace var.
pch=c(19, 17), ### Symbols for levels of trace var.
fixed=TRUE, ### Order by factor order in data
leg.bty = "o") ### Legend box
# Produce an interaction plot for Time*Temperature
interaction.plot(x.factor = Time,
trace.factor = Temperature,
response = dat16$Yield,
fun = mean,
type="b", ### Show plot points as symbols
col=c("forestgreen","orange"), ### Colors for levels of trace var.
pch=c(19, 17), ### Symbols for levels of trace var.
fixed=TRUE, ### Order by factor order in data
leg.bty = "o") ### Legend box