1 Problem

In a process development study on yield, four factors were studied, each at two levels: time (A), concentration (B), pressure (C), and temperature (D).

  1. Display the halfnormal plot for this data and determine which factors appear to be significant.
#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.

  1. 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.
# 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.

2 Complete R Code

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)