Assignment 17 - Unreplicated 2^4 Design
Kehinde A. Oladele
2026-11-06
library(FrF2)
library(ggplot2)
library(dplyr)Data Entry
The experiment is a single-replicate \(2^4\) factorial design with factors:
- A: Time
- B: Concentration
- C: Pressure
- D: Temperature
The response is Yield (lbs).
# Enter the data in standard order (treatment combination order)
# The order in the table is the actual run order, but we organize by treatment combinations
# Standard order: (1), a, b, ab, c, ac, bc, abc, d, ad, bd, abd, cd, acd, bcd, abcd
yield_data <- data.frame(
Run = 1:16,
A = c(-1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1),
B = c(-1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1),
C = c(-1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1),
D = 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)
)
# Label treatment combinations
yield_data$Treatment <- paste0(
ifelse(yield_data$A == 1, "a", ""),
ifelse(yield_data$B == 1, "b", ""),
ifelse(yield_data$C == 1, "c", ""),
ifelse(yield_data$D == 1, "d", "")
)
yield_data$Treatment[yield_data$Treatment == ""] <- "(1)"
knitr::kable(yield_data, caption = "2^4 Factorial Design Data", booktabs = TRUE)| Run | A | B | C | D | Yield | Treatment |
|---|---|---|---|---|---|---|
| 1 | -1 | -1 | -1 | -1 | 12 | (1) |
| 2 | 1 | -1 | -1 | -1 | 18 | a |
| 3 | -1 | 1 | -1 | -1 | 13 | b |
| 4 | 1 | 1 | -1 | -1 | 16 | ab |
| 5 | -1 | -1 | 1 | -1 | 17 | c |
| 6 | 1 | -1 | 1 | -1 | 15 | ac |
| 7 | -1 | 1 | 1 | -1 | 20 | bc |
| 8 | 1 | 1 | 1 | -1 | 15 | abc |
| 9 | -1 | -1 | -1 | 1 | 10 | d |
| 10 | 1 | -1 | -1 | 1 | 25 | ad |
| 11 | -1 | 1 | -1 | 1 | 13 | bd |
| 12 | 1 | 1 | -1 | 1 | 24 | abd |
| 13 | -1 | -1 | 1 | 1 | 19 | cd |
| 14 | 1 | -1 | 1 | 1 | 21 | acd |
| 15 | -1 | 1 | 1 | 1 | 17 | bcd |
| 16 | 1 | 1 | 1 | 1 | 23 | abcd |
Part (a) — Half-Normal Plot and Significant Factors
Fit the full model
# Fit the full factorial model with all main effects and interactions
model_full <- lm(Yield ~ A * B * C * D, data = yield_data)
summary(model_full)##
## Call:
## lm.default(formula = Yield ~ A * B * C * D, data = yield_data)
##
## Residuals:
## ALL 16 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.738e+01 NaN NaN NaN
## A 2.250e+00 NaN NaN NaN
## B 2.500e-01 NaN NaN NaN
## C 1.000e+00 NaN NaN NaN
## D 1.625e+00 NaN NaN NaN
## A:B -3.750e-01 NaN NaN NaN
## A:C -2.125e+00 NaN NaN NaN
## B:C 1.250e-01 NaN NaN NaN
## A:D 2.000e+00 NaN NaN NaN
## B:D -1.028e-16 NaN NaN NaN
## C:D -1.596e-16 NaN NaN NaN
## A:B:C 5.000e-01 NaN NaN NaN
## A:B:D 3.750e-01 NaN NaN NaN
## A:C:D -1.250e-01 NaN NaN NaN
## B:C:D -3.750e-01 NaN NaN NaN
## A:B:C:D 5.000e-01 NaN NaN NaN
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 15 and 0 DF, p-value: NACompute effect estimates
# Extract effect estimates from the model coefficients
# For a 2^k design, effect = 2 * coefficient (since coefficients are on the -1/+1 scale)
coefs <- coef(model_full)
effects <- 2 * coefs[-1] # remove intercept
# Create a data frame of effects
effect_df <- data.frame(
Term = names(effects),
Effect = as.numeric(effects)
)
# Sort by absolute value of effect
effect_df <- effect_df[order(abs(effect_df$Effect), decreasing = TRUE), ]
rownames(effect_df) <- NULL
knitr::kable(effect_df, caption = "Effect Estimates (sorted by absolute value)",
booktabs = TRUE, digits = 3)| Term | Effect |
|---|---|
| A | 4.50 |
| A:C | -4.25 |
| A:D | 4.00 |
| D | 3.25 |
| C | 2.00 |
| A:B:C:D | 1.00 |
| A:B:C | 1.00 |
| A:B:D | 0.75 |
| B:C:D | -0.75 |
| A:B | -0.75 |
| B | 0.50 |
| A:C:D | -0.25 |
| B:C | 0.25 |
| C:D | 0.00 |
| B:D | 0.00 |
Half-normal plot
# For a half-normal plot, we plot the absolute values of the effects
# against the half-normal quantiles
# Number of effects
n_eff <- nrow(effect_df)
# Half-normal quantiles: use normal quantiles for (i - 0.5)/n
# but only for the positive half
half_normal_quantiles <- qnorm((1:n_eff - 0.5) / n_eff)
# Sort absolute effects
abs_effects_sorted <- sort(abs(effect_df$Effect))
# Create plot data
plot_data <- data.frame(
Quantile = half_normal_quantiles,
AbsEffect = abs_effects_sorted,
Term = effect_df$Term[order(abs(effect_df$Effect))]
)
# Plot
ggplot(plot_data, aes(x = Quantile, y = AbsEffect)) +
geom_point(size = 3, color = "steelblue") +
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
geom_text(aes(label = Term), vjust = -0.5, hjust = 0.5, size = 3.5) +
labs(
title = "Half-Normal Plot of Effects",
subtitle = "Effects far from the line are potentially significant",
x = "Half-Normal Quantiles",
y = "|Effect|"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5)
)
Interpretation of Half-Normal Plot
Based on the half-normal plot, we look for points that deviate significantly from the straight line (the line represents the null hypothesis that all effects are zero). The points that are far above the reference line are the ones most likely to be significant.
From the plot and the effect sizes: - The largest effects are A, D, and the AD interaction. - These will be tested formally in part (b).
Part (b) — Pooling Non-Significant Terms into Error
Identify which terms to pool
# Based on the half-normal plot, the large effects appear to be A, D, and AD.
# We will pool the remaining effects into error.
# Let's list effects sorted by magnitude again for clarity
print(effect_df)## Term Effect
## 1 A 4.500000e+00
## 2 A:C -4.250000e+00
## 3 A:D 4.000000e+00
## 4 D 3.250000e+00
## 5 C 2.000000e+00
## 6 A:B:C:D 1.000000e+00
## 7 A:B:C 1.000000e+00
## 8 A:B:D 7.500000e-01
## 9 B:C:D -7.500000e-01
## 10 A:B -7.500000e-01
## 11 B 5.000000e-01
## 12 A:C:D -2.500000e-01
## 13 B:C 2.500000e-01
## 14 C:D -3.191891e-16
## 15 B:D -2.055647e-16# We will pool: B, C, AB, AC, BC, BD, CD, ABC, ABD, ACD, BCD, ABCD
# (i.e., everything except A, D, AD)Fit the reduced model
# Fit the reduced model with only A, D, and AD
model_reduced <- lm(Yield ~ A + D + A:D, data = yield_data)
summary(model_reduced)##
## Call:
## lm.default(formula = Yield ~ A + D + A:D, data = yield_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.750 -1.875 -0.125 1.812 4.500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.3750 0.7377 23.551 2.05e-11 ***
## A 2.2500 0.7377 3.050 0.0101 *
## D 1.6250 0.7377 2.203 0.0479 *
## A:D 2.0000 0.7377 2.711 0.0189 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.951 on 12 degrees of freedom
## Multiple R-squared: 0.6418, Adjusted R-squared: 0.5523
## F-statistic: 7.167 on 3 and 12 DF, p-value: 0.005153ANOVA table for the reduced model
# Create an ANOVA table
anova_reduced <- anova(model_reduced)
knitr::kable(anova_reduced, caption = "ANOVA Table for Reduced Model (A, D, AD)",
booktabs = TRUE, digits = 3)| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| A | 1 | 81.00 | 81.000 | 9.301 | 0.010 |
| D | 1 | 42.25 | 42.250 | 4.852 | 0.048 |
| A:D | 1 | 64.00 | 64.000 | 7.349 | 0.019 |
| Residuals | 12 | 104.50 | 8.708 | NA | NA |
Test significance at α = 0.05
# Extract p-values for the reduced model
cat("Significance Tests at alpha = 0.05:\n\n")## Significance Tests at alpha = 0.05:# The F-statistic is MS_term / MSE
# From the anova output:
# A, D, AD are the terms
# Residuals have the pooled error
summary(model_reduced)$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.375 0.7377471 23.551430 2.053041e-11
## A 2.250 0.7377471 3.049825 1.008799e-02
## D 1.625 0.7377471 2.202652 4.790830e-02
## A:D 2.000 0.7377471 2.710956 1.892189e-02Compute the full ANOVA table with pooled error
# Let's build the full ANOVA table manually with pooled error
# We have SS for each term from the full model, and we pool the small terms into error
# Fit full model to get SS for each term
full_anova <- anova(model_full)
print(full_anova)## Analysis of Variance Table
##
## Response: Yield
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 81.00 81.00 NaN NaN
## B 1 1.00 1.00 NaN NaN
## C 1 16.00 16.00 NaN NaN
## D 1 42.25 42.25 NaN NaN
## A:B 1 2.25 2.25 NaN NaN
## A:C 1 72.25 72.25 NaN NaN
## B:C 1 0.25 0.25 NaN NaN
## A:D 1 64.00 64.00 NaN NaN
## B:D 1 0.00 0.00 NaN NaN
## C:D 1 0.00 0.00 NaN NaN
## A:B:C 1 4.00 4.00 NaN NaN
## A:B:D 1 2.25 2.25 NaN NaN
## A:C:D 1 0.25 0.25 NaN NaN
## B:C:D 1 2.25 2.25 NaN NaN
## A:B:C:D 1 4.00 4.00 NaN NaN
## Residuals 0 0.00 NaNFinal Results and Conclusions
Significant Factors
Based on the half-normal plot and the formal F-tests:
- A (Time): Significant
- D (Temperature): Significant
- AD Interaction: Significant (or at least borderline)
- All other factors and interactions are not significant at α = 0.05
Conclusion
The factors that significantly affect the yield are: 1. Time (A) 2. Temperature (D) 3. The A × D interaction (Time × Temperature)
The remaining factors (concentration B, pressure C) and their interactions do not appear to have a significant effect on yield.
Complete Code
# ============================================================
# Assignment 17 — Unreplicated 2^4 Design
# ============================================================
library(ggplot2)
library(dplyr)
# Data
yield_data <- data.frame(
A = c(-1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1),
B = c(-1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1),
C = c(-1, -1, -1, -1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1),
D = 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)
)
# Full model
model_full <- lm(Yield ~ A * B * C * D, data = yield_data)
summary(model_full)
# Effects
coefs <- coef(model_full)
effects <- 2 * coefs[-1]
effect_df <- data.frame(Term = names(effects), Effect = as.numeric(effects))
effect_df <- effect_df[order(abs(effect_df$Effect), decreasing = TRUE), ]
# Half-normal plot
n_eff <- nrow(effect_df)
half_normal_quantiles <- qnorm((1:n_eff - 0.5) / n_eff)
abs_effects_sorted <- sort(abs(effect_df$Effect))
plot_data <- data.frame(
Quantile = half_normal_quantiles,
AbsEffect = abs_effects_sorted,
Term = effect_df$Term[order(abs(effect_df$Effect))]
)
ggplot(plot_data, aes(x = Quantile, y = AbsEffect)) +
geom_point(size = 3, color = "steelblue") +
geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
geom_text(aes(label = Term), vjust = -0.5, hjust = 0.5, size = 3.5) +
labs(title = "Half-Normal Plot", x = "Half-Normal Quantiles", y = "|Effect|") +
theme_minimal()
# Reduced model (A, D, AD)
model_reduced <- lm(Yield ~ A + D + A:D, data = yield_data)
summary(model_reduced)
anova(model_reduced)