實驗設計與反應曲面法 Design of Experiments & Response Surface Methodology — with Base R

研發會議上主管說：「下個月要交出最佳萃取條件，實驗預算只夠 20 次。」OFAT 走不完、瞎試沒根據——這一章就是你的預算內攻略。

14.0 DOE 入門：先認識「實驗設計」這件事

本章所有內容建立在一個前提上：實驗的「設計」比「執行」更早決定成敗。在學統計語言之前，先用一個例子把觀念建立起來。

有了這些觀念，下面正式比較 OFAT 與 DOE 的差別。

14.1 為什麼是 DOE：OFAT 的陷阱

14.2 2³ 全因子 + 中心點：效應與曲率

8 個因子點（±1 組合全跑）＋ 4 個中心點（0,0,0 重複）。設計矩陣與模擬結果：

配一階＋交互作用模型 lm(y ~ (A+B+C)^2)：

圖 14-1 左：A×B 交互作用（兩線不平行）與中心點曲率（綠色菱形遠高於兩線）；右：CCD 二階模型的等高線與駐點

曲率檢定：中心點的殺手鐧

2 水準設計只能配「直線」——但中心點實測平均 60.2%，因子點平均只有 43.0%， 差了 17 個百分點。用中心點的純誤差 s = 0.86 做檢定：

14.3 RSM 入門：把「找最大值」變成「找山頂」

14.2 的曲率檢定告訴我們「曲面是彎的」——這正是 RSM（Response Surface Methodology，反應曲面法）登場的時機。

14.3.1 中心複合設計 (CCD)：二階模型與駐點

CCD 在 2³ 因子設計上補軸點（每因子 ±α，其他因子歸零）與更多中心點：

二階模型 lm(y ~ (A+B+C)^2 + I(A^2)+I(B^2)+I(C^2)) 的配適結果：

三個平方項全為負 → 開口朝下 → 存在極大值（呼應曲率檢定）。把模型寫成矩陣形式 ŷ = b₀ + x′b + x′Bx（B 對角線放平方項、非對角線放交互係數的一半），駐點有閉式解：

14.3.2 模型檢驗：ANOVA 表與失擬檢定

把駐點交給操作員之前，先問模型兩個問題：① 哪些項真的在出力？（ANOVA 表，Ch16 的「把總變異拆開」） ② 二階模型的「形狀」夠不夠？（失擬檢定，Ch17 17.8–17.10）。

① anova(fit2)：一階、交互作用、二次項各吃掉多少 SS

20 個 run 的總變異 SS_總 = 1894.66（df = 19）。anova(fit2) 把它逐項拆開（完整輸出見本章末的 R Console），依項目類型加總：

• 每一列的讀法和 Ch16 一樣：F = 該項的 MS ÷ MS_殘差；H₀ 是「該項係數＝0」。例如 A：348.66 / 1.18 = 295.3、p = 9.4×10⁻⁹。 每項只有 1 個 df，所以這個 F 就是係數 t 值的平方（A 的 t = 17.19，17.19² ≈ 295）。
• A:C、B:C 不顯著——和 14.2 的「可忽略」一致（教學模擬的真實模型裡本來就沒有這兩項）。
• 注意順序：anova() 給的是「依序放入」的 SS（R 會把 I(A^2) 等平方項排在交互作用之前）。 CCD 的一階項與交互作用項彼此正交，SS 不受順序影響；但三個平方項彼此不完全正交，個別 SS 會隨順序略變， 三項合計 1318.12 則不變——所以二次項要「整組」讀。

② 失擬檢定：6 個中心點的真正用途

20 個 run 裡只有 15 個不同的設計點；唯一有重複的就是 6 個一模一樣的中心點。它們不是浪費—— 同一條件重做 6 次的差異與模型無關，就是純誤差（s = 1.054，SS = 5 × 1.054² = 5.55，df = 6 − 1 = 5）。 有了它，就能照 Ch17 的做法把殘差 SS 拆成兩塊：

同樣的檢定也對「只有一階＋交互作用」的模型做一次，對照就很清楚：

兩個模型的殘差 SS 相差 1329.93 − 11.81 = 1318.12，正好就是上面 ANOVA 表的「二次項」那一塊——一階模型的失擬幾乎全是沒配到的曲率。 這也是 14.2「曲率檢定」的另一種問法：那裡用 4 個中心點比平均，這裡用 6 個中心點比 SS，結論相同。

14.4 預測區間與驗證實驗

報告最適條件時，單點預測必須附預測區間（呼應 Ch2 的 CI 與 Ch8 的報告規範）：

14.5 套件選讀：CRAN 的 DoE 生態系

前四節證明了：因子設計與 RSM 的原理完全在 base R 之內。實務上則用套件省下建矩陣與手算的時間 （先 install.packages(c("rsm", "FrF2", "desirability"))，再把腳本的 HAVE_PKGS 改 TRUE）：

# ====== Ch14 DOE & RSM（完整版請下載 ch14_doe_rsm.R）======
# 情境：超音波輔助萃取茶葉多酚，DPPH(%) ← A乙醇濃度 / B溫度 / C時間
# --- Part 1. 2^3 全因子 + 中心點：效應與曲率（base R）---
true_y <- function(A, B, C)                 # 教學模擬的「真實模型」
  60 + 5*A + 3*B + 2*C - 7*A^2 - 5*B^2 - 6*C^2 + 2.5*A*B
fac <- expand.grid(A = c(-1,1), B = c(-1,1), C = c(-1,1))
set.seed(14)
fac$y <- round(true_y(fac$A, fac$B, fac$C) + rnorm(8, 0, 1.2), 1)
ctr <- data.frame(A=0, B=0, C=0, y=round(rnorm(4, 60, 1.2), 1))  # 中心點x4
dat1 <- rbind(fac, ctr)
fit1 <- lm(y ~ (A + B + C)^2, dat1)
round(coef(fit1), 2)
#> A=5.64 B=3.34 C=1.64 A:B=2.04（效應 = 2*係數；A:C、B:C 可忽略）
# 曲率檢定：中心點平均 vs 因子點平均，純誤差來自中心點重複
t_curv <- (mean(ctr$y) - mean(fac$y)) / (sd(ctr$y) * sqrt(1/4 + 1/8))
round(c(mean_fac = mean(fac$y), mean_ctr = mean(ctr$y),
        t = t_curv, p = 2 * pt(-abs(t_curv), 3)), 4)
#> 43.04 vs 60.18, t = 32.5, p = 0.0001 -> 曲面是彎的！需要 RSM

# --- Part 2. 中心複合設計 CCD：二階模型 + 駐點 ---
alpha <- 2^(3/4)                            #> 1.682（k=3 旋轉性設計）
axi <- data.frame(A = c(-alpha, alpha, rep(0, 4)),
                  B = c(0, 0, -alpha, alpha, 0, 0),
                  C = c(0, 0, 0, 0, -alpha, alpha))
ccd <- rbind(fac[, c("A","B","C")], axi,
             data.frame(A=0, B=0, C=0)[rep(1, 6), ])   # 8+6+6 = 20 runs
set.seed(15)
ccd$y <- round(true_y(ccd$A, ccd$B, ccd$C) + rnorm(20, 0, 1.2), 1)
fit2 <- lm(y ~ (A + B + C)^2 + I(A^2) + I(B^2) + I(C^2), ccd)
round(coef(fit2), 2)
#> A^2=-7.08 B^2=-4.71 C^2=-6.04 全負 -> 開口朝下，有極大值
# 駐點 x0 = -0.5 * B^-1 * b（B 對角=平方項，非對角=交互係數/2）
b <- coef(fit2); blin <- b[c("A", "B", "C")]
Bmat <- matrix(c(b["I(A^2)"], b["A:B"]/2, b["A:C"]/2,
                 b["A:B"]/2, b["I(B^2)"], b["B:C"]/2,
                 b["A:C"]/2, b["B:C"]/2, b["I(C^2)"]), 3, byrow = TRUE)
x0 <- as.numeric(-0.5 * solve(Bmat, blin))
eigen(Bmat)$values                          #> -7.82 -5.98 -4.02 全負=極大值
y0 <- as.numeric(b[1] + t(x0) %*% blin + t(x0) %*% Bmat %*% x0)
round(c(A = x0[1], B = x0[2], C = x0[3], DPPH = y0), 2)
#> (0.44, 0.45, 0.13)，預測 61.8%
round(c(乙醇 = 50 + 10*x0[1], 溫度 = 50 + 10*x0[2],
        時間 = 22.5 + 7.5*x0[3]), 1)
#> 乙醇 54.4% / 溫度 54.5°C / 超音波 23.5 min

# --- Part 2b. 模型檢驗：ANOVA 表 + 失擬檢定（呼應 Ch16/Ch17）---
aov2 <- anova(fit2)
aov2                                        # 每列 H0：該項係數 = 0
ss <- setNames(aov2[["Sum Sq"]], rownames(aov2))
round(c(一階 = sum(ss[c("A", "B", "C")]),
        二次 = sum(ss[c("I(A^2)", "I(B^2)", "I(C^2)")]),
        交互 = sum(ss[c("A:B", "A:C", "B:C")]),
        殘差 = unname(ss["Residuals"]), 總和 = sum(ss)), 2)
#> 一階 494.40 / 二次 1318.12 / 交互 70.33 / 殘差 11.81 / 總和 1894.66
# 失擬檢定：只有 6 個中心點有重複 -> 純誤差全靠它們 (df = 5)
ccd$pt <- factor(paste(ccd$A, ccd$B, ccd$C))   # 每個設計點一個水準：15 個
fit_pt <- lm(y ~ pt, data = ccd)               # 每點一個平均 -> 殘差 = 純誤差
fit1c  <- lm(y ~ (A + B + C)^2, data = ccd)    # 對照：只有一階+交互作用
round(sd(ccd$y[15:20]), 3)                  #> 1.054 = 中心點的 s（純誤差）
anova(fit1c, fit_pt)   #> F = 149.05, p = 1.6e-05 -> 一階模型失擬顯著（有曲率）
anova(fit2, fit_pt)    #> F = 1.13, p = 0.45 -> 二階模型失擬不顯著，夠用

# --- Part 3. 預測區間與驗證實驗 ---
pred <- predict(fit2, data.frame(A = x0[1], B = x0[2], C = x0[3]),
                interval = "prediction")
round(pred, 1)                              #> 61.8 (59.2 ~ 64.4) 95% PI
set.seed(16)
verify <- round(true_y(x0[1], x0[2], x0[3]) + rnorm(3, 0, 1.2), 1)
#> 62.4 61.7 63.2 -> 平均落在區間內 = 驗證通過（呼應 Ch2/Ch8）

# --- Part 4. 套件選讀（先 install.packages，再改 HAVE_PKGS <- TRUE）---
# HAVE_PKGS <- FALSE
# if (HAVE_PKGS) {
#   library(rsm)                             # Lenth (2009) JSS 32:7
#   fit_rsm <- rsm(y ~ SO(A, B, C), data = ccd)   # 二階配適
#   summary(fit_rsm)$canonical               # 駐點 + canonical 分析
#   steepest(fit_rsm)                        # 最陡上升路徑表
#   contour(fit_rsm, ~ A + B, image = TRUE)  # 一行等高線+填色
#   library(FrF2)                            # Grömping (2014) JSS 56:1
#   FrF2(8, 4)                               # 2^(4-1) 半因子：4因子只用8 run
#   pb(12, 5)                                # Plackett-Burman 篩 5 變因
#   library(desirability)                    # 多反應最佳化
#   d1 <- dMax(40, 65)                       # DPPH 越高越好
#   d2 <- dMin(30, 60)                       # 溶劑用量越少越好
# }

(Intercept)           A           B           C         A:B         A:C 
      48.75        5.64        3.34        1.64        2.04        0.09 
        B:C 
      -0.41 
mean_fac mean_ctr        t        p 
 43.0375  60.1750  32.4776   0.0001 
(Intercept)           A           B           C      I(A^2)      I(B^2) 
      59.94        5.05        2.86        1.58       -7.08       -4.71 
     I(C^2)         A:B         A:C         B:C 
      -6.04        2.88       -0.52        0.50 
[1] -4.024439 -5.975244 -7.822544
    A     B     C  DPPH 
 0.44  0.45  0.13 61.80 
乙醇 溫度 時間 
54.4 54.5 23.5 
Analysis of Variance Table

Response: y
          Df Sum Sq Mean Sq  F value    Pr(>F)    
A          1 348.66  348.66 295.3367 9.402e-09 ***
B          1 111.55  111.55  94.4930 2.059e-06 ***
C          1  34.19   34.19  28.9624 0.0003094 ***
I(A^2)     1 547.39  547.39 463.6737 1.041e-09 ***
I(B^2)     1 245.86  245.86 208.2561 5.067e-08 ***
I(C^2)     1 524.88  524.88 444.6092 1.279e-09 ***
A:B        1  66.13   66.13  56.0123 2.101e-05 ***
A:C        1   2.20    2.20   1.8678 0.2016677    
B:C        1   2.00    2.00   1.6941 0.2222425    
Residuals 10  11.81    1.18                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   一階    二次    交互    殘差    總和 
 494.40 1318.12   70.33   11.81 1894.66 
[1] 1.054
Analysis of Variance Table

Model 1: y ~ (A + B + C)^2
Model 2: y ~ pt
  Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
1     13 1329.93                                  
2      5    5.55  8    1324.4 149.05 1.612e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Analysis of Variance Table

Model 1: y ~ (A + B + C)^2 + I(A^2) + I(B^2) + I(C^2)
Model 2: y ~ pt
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1     10 11.8054                           
2      5  5.5533  5    6.2521 1.1258 0.4498
   fit  lwr  upr
1 61.8 59.2 64.4

📊 Excel 也會算：DOE 效應與 RSM

14.6 參考文獻

1. Lenth, R. V. Response-Surface Methods in R, Using the rsm Package. J. Stat. Software 2009, 32 (7). DOI: 10.18637/jss.v032.i07
2. Grömping, U. R Package FrF2 for Creating and Analyzing Fractional Factorial 2-Level Designs. J. Stat. Software 2014, 56 (1). DOI: 10.18637/jss.v056.i01
3. Tanaka, E.; Amaliah, D. Current State and Prospects of R-Packages for the Design of Experiments. 2022. arXiv:2206.07532（CRAN ExperimentalDesign 任務視圖 105 套件的量化總體檢）.
4. Myers, R. H.; Montgomery, D. C.; Anderson-Cook, C. M. Response Surface Methodology: Process and Product Optimization Using Designed Experiments, 4th ed.; Wiley, 2016.（RSM 的聖經：CCD/BBD、canonical 分析、多反應最佳化）.
5. Smith, J. S. Evaluation of Analytical Data. In Nielsen’s Food Analysis, 6th ed.; Ismail, B. P., Nielsen, S. S., Eds.; Springer, 2024; Chapter 4.（本課程數據品管統計之主軸教材）.