1 Problem

A process development study investigated yield using four two-level factors:

The response is yield in pounds. The experiment is an unreplicated full \(2^4\) factorial design. The table gives the standard treatment combination (Run Number), the randomized actual run order, factor settings, and observed yield.

2 1. Enter and inspect the data

# Responses are entered in standard factorial Run Number order.
# Actual Run Order is recorded separately and does not change the
# treatment-to-response pairing shown in the assignment table.

actual_order <- c(5, 9, 8, 13, 3, 7, 14, 1,
                  6, 11, 2, 15, 4, 16, 10, 12)

yield <- c(12, 18, 13, 16, 17, 15, 20, 15,
           10, 25, 13, 24, 19, 21, 17, 23)

# Standard 2^4 order: A changes fastest, then B, C, and D.
design <- expand.grid(
  A = c(-1, 1),
  B = c(-1, 1),
  C = c(-1, 1),
  D = c(-1, 1)
)

data <- data.frame(
  Run_Number = 1:16,
  Actual_Run_Order = actual_order,
  design,
  Yield_lbs = yield
)

knitr::kable(data, caption = "Unreplicated 2^4 factorial experiment")
Unreplicated 2^4 factorial experiment
Run_Number Actual_Run_Order A B C D Yield_lbs
1 5 -1 -1 -1 -1 12
2 9 1 -1 -1 -1 18
3 8 -1 1 -1 -1 13
4 13 1 1 -1 -1 16
5 3 -1 -1 1 -1 17
6 7 1 -1 1 -1 15
7 14 -1 1 1 -1 20
8 1 1 1 1 -1 15
9 6 -1 -1 -1 1 10
10 11 1 -1 -1 1 25
11 2 -1 1 -1 1 13
12 15 1 1 -1 1 24
13 4 -1 -1 1 1 19
14 16 1 -1 1 1 21
15 10 -1 1 1 1 17
16 12 1 1 1 1 23 
cat("Number of runs:", nrow(data), "\n")
## Number of runs: 16
cat("Grand mean yield:", mean(data$Yield_lbs), "lbs\n")
## Grand mean yield: 17.375 lbs

3 2. Calculate all factorial effects

For a \(2^4\) design, each effect estimate is its factorial contrast divided by \(2^{4-1}=8\). The 15 effects include four main effects, six two-factor interactions, four three-factor interactions, and one four-factor interaction.

factor_names <- c("A", "B", "C", "D")
effect_names <- c(
  "A", "B", "C", "D",
  "AB", "AC", "AD", "BC", "BD", "CD",
  "ABC", "ABD", "ACD", "BCD", "ABCD"
)

# Build each interaction contrast from the coded factor columns.
contrast_sign <- function(term, x) {
  term_factors <- strsplit(term, split = "")[[1]]
  Reduce(`*`, lapply(x[term_factors], as.numeric))
}

effect_estimates <- vapply(
  effect_names,
  function(term) {
    sum(contrast_sign(term, design) * data$Yield_lbs) / 8
  },
  numeric(1)
)

# For a 2^4 design, the sum of squares for each one-df effect is
# 2^(4-2) * effect^2 = 4 * effect^2.
effect_table <- data.frame(
  Effect = names(effect_estimates),
  Estimate = as.numeric(effect_estimates),
  Absolute_Effect = abs(as.numeric(effect_estimates)),
  SS = 4 * as.numeric(effect_estimates)^2,
  row.names = NULL
)

effect_table <- effect_table[order(-effect_table$Absolute_Effect), ]
knitr::kable(effect_table, digits = 4,
             caption = "Estimated factorial effects, ordered by absolute magnitude")
Estimated factorial effects, ordered by absolute magnitude
Effect Estimate Absolute_Effect SS
1 A 4.50 4.50 81.00
6 AC -4.25 4.25 72.25
7 AD 4.00 4.00 64.00
4 D 3.25 3.25 42.25
3 C 2.00 2.00 16.00
11 ABC 1.00 1.00 4.00
15 ABCD 1.00 1.00 4.00
5 AB -0.75 0.75 2.25
12 ABD 0.75 0.75 2.25
14 BCD -0.75 0.75 2.25
2 B 0.50 0.50 1.00
8 BC 0.25 0.25 0.25
13 ACD -0.25 0.25 0.25
9 BD 0.00 0.00 0.00
10 CD 0.00 0.00 0.00

4 3. Half-normal plot and effect screening

A half-normal plot is used to screen the effects because the experiment is unreplicated and has no direct pure-error estimate. Effects that follow the small-effect reference pattern are treated as inactive; effects that stand apart from that pattern are considered potentially active.

# Order absolute effects from smallest to largest.
ordered_effects <- effect_table[order(effect_table$Absolute_Effect), ]
n_effects <- nrow(ordered_effects)

# Half-normal plotting scores.
plotting_positions <- (seq_len(n_effects) - 0.5) / n_effects
half_normal_scores <- qnorm((1 + plotting_positions) / 2)

# Screening based on the clear separation in the half-normal plot.
# Five effects depart from the near-linear cluster of smaller effects.
# The visually separated effects are A, C, D, AC, and AD.
screened_active <- c("A", "C", "D", "AC", "AD")

plot(
  half_normal_scores,
  ordered_effects$Absolute_Effect,
  pch = 19,
  xlab = "Half-normal score",
  ylab = "Absolute effect estimate",
  main = "Half-normal Plot of the 15 Factorial Effects"
)
text(
  half_normal_scores,
  ordered_effects$Absolute_Effect,
  labels = ordered_effects$Effect,
  pos = 4,
  cex = 0.85
)

# Add a reference line fitted to the 11 effects screened as inactive.
inactive_for_line <- ordered_effects$Effect %in% setdiff(effect_names, screened_active)
inactive_abs <- ordered_effects$Absolute_Effect[inactive_for_line]
inactive_scores <- half_normal_scores[inactive_for_line]
reference_fit <- lm(inactive_abs ~ inactive_scores)

# Draw the reference only over the score range occupied by inactive effects.
ref_range <- range(inactive_scores)
ref_x <- seq(ref_range[1], ref_range[2], length.out = 100)
ref_y <- pmax(0, as.numeric(coef(reference_fit)[1] +
                             coef(reference_fit)[2] * ref_x))
lines(ref_x, ref_y, lty = 2, lwd = 2)

cat("Effects identified as potentially active from the half-normal plot:\n")
## Effects identified as potentially active from the half-normal plot:
print(screened_active)
## [1] "A"  "C"  "D"  "AC" "AD"

Based on the half-normal plot, the effects that stand apart from the small-effect pattern are A, C, D, AC, and AD. The remaining 10 effects are pooled to estimate experimental error in the next step.

5 4. Pool the effects screened as inactive and test the remaining effects

The 10 effects screened as inactive are pooled into error. Each effect has one degree of freedom, so the pooled error has 10 degrees of freedom. The five retained effects are then tested using this pooled mean square error.

inactive_effects <- setdiff(effect_names, screened_active)
active_effects <- screened_active

# Sum of squares for active and inactive effects.
ss_by_effect <- setNames(effect_table$SS, effect_table$Effect)

SSE <- sum(ss_by_effect[inactive_effects])
df_error <- length(inactive_effects)
MSE <- SSE / df_error

active_anova <- data.frame(
  Source = active_effects,
  DF = 1,
  SS = as.numeric(ss_by_effect[active_effects]),
  MS = as.numeric(ss_by_effect[active_effects])
)

active_anova$F_value <- active_anova$MS / MSE
active_anova$P_value <- pf(
  active_anova$F_value,
  df1 = 1,
  df2 = df_error,
  lower.tail = FALSE
)
active_anova$Decision <- ifelse(
  active_anova$P_value < alpha,
  "Significant",
  "Not significant"
)

pooled_error_row <- data.frame(
  Source = "Pooled error (10 inactive effects)",
  DF = df_error,
  SS = SSE,
  MS = MSE,
  F_value = NA_real_,
  P_value = NA_real_,
  Decision = "Error estimate"
)

anova_table <- rbind(active_anova, pooled_error_row)

knitr::kable(anova_table, digits = 6,
             caption = "ANOVA after pooling the 11 effects screened as inactive")
ANOVA after pooling the 11 effects screened as inactive
Source DF SS MS F_value P_value Decision
A 1 81.00 81.000 49.846154 0.000035 Significant
C 1 16.00 16.000 9.846154 0.010549 Significant
D 1 42.25 42.250 26.000000 0.000465 Significant
AC 1 72.25 72.250 44.461538 0.000056 Significant
AD 1 64.00 64.000 39.384615 0.000092 Significant
Pooled error (10 inactive effects) 10 16.25 1.625 NA NA Error estimate 
cat("Pooled error degrees of freedom:", df_error, "\n")
## Pooled error degrees of freedom: 10
cat("Pooled error sum of squares:", round(SSE, 6), "\n")
## Pooled error sum of squares: 16.25
cat("Pooled error mean square:", round(MSE, 6), "\n")
## Pooled error mean square: 1.625
# Show the total corrected sum of squares as a check.
SST <- sum((data$Yield_lbs - mean(data$Yield_lbs))^2)
SS_effects <- sum(effect_table$SS)
cat("Total corrected SS:", round(SST, 6), "\n")
## Total corrected SS: 291.75
cat("Sum of all 15 factorial effect SS:", round(SS_effects, 6), "\n")
## Sum of all 15 factorial effect SS: 291.75
stopifnot(abs(SST - SS_effects) < 1e-8)

6 5. Conclusions

At the 0.05 significance level, the pooled-error tests identify the following effects as statistically significant:

The other 10 effects were screened as inactive using the half-normal plot and pooled to estimate error. The retained effects have p-values below 0.05 when tested against the pooled error mean square.

The positive C effect indicates that increasing pressure from its low level to its high level increases mean yield, averaged over the other factors. The negative AC effect (estimate = -4.25 lbs) indicates that the combined effect of time and pressure is opposite in direction to what would be expected from simply adding their separate main effects. The positive A, C, D, and AD estimates are 4.50, 2.00, 3.25, and 4.00 lbs, respectively. The positive AD effect indicates a positive time-by-temperature interaction. Since interactions AC and AD are significant, the effects of time should be interpreted in relation to pressure and temperature rather than only as isolated main effects.

Overall conclusion: The half-normal plot indicates that A, C, D, AC, and AD are the potentially active effects. After pooling the other effects, all five retained effects are statistically significant at \(\alpha=0.05\).

7 Reproducibility check

stopifnot(nrow(data) == 16)
stopifnot(length(effect_estimates) == 15)
stopifnot(df_error == 10)
stopifnot(all(active_anova$P_value < alpha))
cat("All data, effect-count, pooled-error, and significance checks passed.\n")
## All data, effect-count, pooled-error, and significance checks passed.