In a process development study on yield, four factors were studied, each at two levels: time (A), concentration (B), pressure (C), and temperature (D).
#install.packages("FrF2")
library("FrF2")
# Data
design_data <- data.frame(
Run = 1:16,
ActualRunOrder = c(5, 9, 8, 13, 3, 7, 14, 1, 6, 11, 2, 15, 4, 16, 10, 12),
A = factor(c(-1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1)),
B = factor(c(-1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1)),
C = factor(c(-1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1)),
D = factor(c(-1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1)),
Yield = c(12, 18, 13, 16, 17, 15, 20, 15, 10, 25, 13, 24, 19, 21, 17, 23)
)
# Display
print(design_data)
## Run 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
# Fit linear model with numeric levels (-1, 1) to extract effects
A_num <- as.numeric(as.character(design_data$A))
B_num <- as.numeric(as.character(design_data$B))
C_num <- as.numeric(as.character(design_data$C))
D_num <- as.numeric(as.character(design_data$D))
full_model <- lm(Yield ~ A_num * B_num * C_num * D_num, data = design_data)
# Extract effects (2 * coefficients, excluding intercept)
effects <- 2 * coef(full_model)[-1]
# Half-normal plot Daniel's method (DanielPlot from FrF2)
DanielPlot(full_model, main = "Half-Normal Plot of Effects")
Based on the half-normal plot, the significant effects are: A, C, D, AC, and AD. All other terms are pooled into error.
# ANOVA with pooled error
anova_model <- lm(Yield ~ A + C + D + A:C + A:D, data = design_data)
anova_table <- anova(anova_model)
print(anova_table)
## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 81.00 81.000 49.8462 3.456e-05 ***
## C 1 16.00 16.000 9.8462 0.0105485 *
## D 1 42.25 42.250 26.0000 0.0004647 ***
## A:C 1 72.25 72.250 44.4615 5.583e-05 ***
## A:D 1 64.00 64.000 39.3846 9.193e-05 ***
## Residuals 10 16.25 1.625
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conclusion: Significant Factors at alpha = 0.05
Factor A (Time): p-value < 0.05, therefore it is statistically significant.
Factor C (Pressure): p-value < 0.05, therefore it is statistically significant.
Factor D (Temperature): p-value < 0.05, therefore it is statistically significant.
The interaction A:C (Time x Pressure): p-value < 0.05, therefore it is statistically significant.
The interaction A:D (Time x Temperature): p-value < 0.05, therefore it is statistically significant.
knitr::opts_chunk$set(echo = TRUE, warning=FALSE, message = FALSE)
#install.packages("FrF2")
library("FrF2")
# Data
design_data <- data.frame(
Run = 1:16,
ActualRunOrder = c(5, 9, 8, 13, 3, 7, 14, 1, 6, 11, 2, 15, 4, 16, 10, 12),
A = factor(c(-1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1)),
B = factor(c(-1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1)),
C = factor(c(-1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1)),
D = factor(c(-1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1)),
Yield = c(12, 18, 13, 16, 17, 15, 20, 15, 10, 25, 13, 24, 19, 21, 17, 23)
)
# Display
print(design_data)
# Fit linear model with numeric levels (-1, 1) to extract effects
A_num <- as.numeric(as.character(design_data$A))
B_num <- as.numeric(as.character(design_data$B))
C_num <- as.numeric(as.character(design_data$C))
D_num <- as.numeric(as.character(design_data$D))
full_model <- lm(Yield ~ A_num * B_num * C_num * D_num, data = design_data)
# Extract effects (2 * coefficients, excluding intercept)
effects <- 2 * coef(full_model)[-1]
# Half-normal plot Daniel's method (DanielPlot from FrF2)
DanielPlot(full_model, main = "Half-Normal Plot of Effects")
# ANOVA with pooled error
anova_model <- lm(Yield ~ A + C + D + A:C + A:D, data = design_data)
anova_table <- anova(anova_model)
print(anova_table)