用R解析DOE實驗數據令人驚豔

以Excel工作表05的2^4 四因子完整因子設計為例說明R的解析過程

1 用Excel數據輸入並儲存：
雖然RCommander的功能表上有[Edit data set]可以直接輸入Data，但用慣Excel的我覺得很不好用，因此還是用Excel輸入然後儲存，檔名與位址 d:\DOE_R\DOE05.csv，下圖是以windows記事本開啟的部分內容，注意首列式標題，紀載前4項為因子，末項Y為實驗結果，用R語言特別注意文字大小寫不能混錯。

(註：以下R指令中要指定名稱或數值的操作符號原本是『 < - 』，但鍵入比較麻煩，最近R的版本已經接受『 = 』符號代替，後文我都使用等號。

2 讓R讀取數據csv檔案，使用2個指定讓R建立數據組

tfile="D:/DOE_R/DOE05.csv" data1= read.csv(tfile, header = TRUE, sep = ",")

用Copy and Paste將2指令複製到RCommander的R Script方框內，因為一次2個指令，需先選取後再按[Submit]鍵入R中，一般未出現紅色警語大致表示R已接受指令，若要確認數據是否順利讀入，可打上指令data1後按[Submit]，就會列出該csv檔的資料內容。(注意指令中『/』為目錄符號，與windows慣用的『\』有差別。)

3 多因子將線性模型進行首次解析，輸入底下指令

lm.out =lm(Y ~ A*B*C*D, data=data1) lm.out

指令lm是進行縣性模型分析，結果放在lm.out物件內，本例為4個因子故以A*B*C*D表示lm迴歸式中的formula，A*B*C*D表示迴歸式的Terms將包含4個主因子、10個二因子交互作用、6個三因子交互作用、1個因子交互作用等項目，這種表示比用JMP、Minitab來得簡潔，我很喜歡。R輸出經Least Square後得到迴歸式係數，截距項為Y值總平均，除截距項外還有15個係數，此係數的+-表示效應方向，係數值的4倍就是效應值，可用指令effects(lm.out)看到Terms各項的效應值。

 居於效應的稀疏性法則，原本希望經檢定找出有顯著性的因子作為關鍵的因子，但迴歸式的Terms包含所有項目，所以誤差項自由度為0，因此ANOVA或母迴歸係數都無法檢定，這情況下傳統的作法是逐步從Terms中剃除係數(效應)較小者作為誤差項，一般高次交互作用是首要目標，這種方式稱為誤差項合併(pooling)，藉誤差項自由度提高而增加檢出顯著性的可能，但相當花費時間。

4 繪製效應常態機率圖(Half-Normal Plot of Effects)
專業的統計軟體常繪製效應常態機率圖，標示著具有顯著性的主效應或交互作用等Terms項目，這樣不必按部就班誤差項合併而可以節省很多時間，R軟體有這個指令稱為DanielPlot，另外BsMD套件有LenthPlot可以較詳細分析潛在顯著性，所以我將2個圖放在一起，但第一次執行LenthPlot指令前需要事先安裝BsMD套件(Package，安裝方法如同安裝Rcmd套件一樣。

library(BsMD) par = par(oma=c(0,2,0,0), mfrow=c(2,1)) DanielPlot(lm.out, half = TRUE, main = "Half-Normal Plot of Effects") LenthPlot(lm.out, main = "Lenth's Plot")

輸出二個圖

註：因為後面主效應與交互作用圖會使用FrF2套件，此套件修改BsMD套件的DanielPlot變得較方便，不顯著的因子名稱也自動隱藏，但BsMD的LenthPlot可以讓我確認PSE值。

library(FrF2) par = par(oma=c(0,2,0,0), mfrow=c(2,1)) DanielPlot(lm.out, alpha = 0.05, half = TRUE, main = "Half-Normal Plot of Effects") DanielPlot(lm.out, alpha = 0.05, half = F, main = "Normal Plot of Effects")

其實效應常態機率圖只要一個就可以，建議選擇Half-Normal Plot of Effects就可以

FrF2套件輸出二個圖

有關輸出圖形的解讀：
1) 基於誤差項自由度為0，無法得到MSe而不能進行因子顯著性檢定，因此統計上出現名為Lenth’s PSE的虛擬誤差項(PSE)的方法可以初步估計各因子的顯著性，PSE的估計是將15個效應值取決對值後計算其中位數(Median)，Lenth取1.5倍中位數作為PSE。
2) DanielPlot(lm.out, half = TRUE, main = "Half-Normal Plot of Effects")是將所有迴歸後的15個效應值之絕對值若大於2.5×PSE就歸為具有顯著性，並以常態機率圖方式作圖，上面Half-Normal Plot of Effect自大而小為A、AC、AD、D、C，應可視為有顯著性，也就是其餘因子可自Terms剃除，然後再次迴歸。
3) 若想區分效應的正負方向，則可將指令改為DanielPlot(lm.out, main = "Normal Plot of Effects")。
4) 下圖Lenth’s Plot係用直條圖表式效應大小，超過ME虛線外表是具有顯著性。

5 剔除不顯著項目再次迴歸分析，輸入底下指令

lm.out =lm(Y ~ A*C+A*D, data=data1) anova(lm.out) summary(lm.out)

指令中的formula Y ~ A*C+A*D也可以這樣寫Y~A+C+D+A:C+A:D，二者是一致的

輸出二個表

1) ANOVA表

表中A:C、A:D表示AC、AD交互作用(R軟體表示法)

2) 母迴歸係數顯著性檢定表

與Excel或其他軟體迴歸分析輸出內容是一致的，共3個部分

其一：迴歸統計，上表的下方可看到R-squared=0.966，Adjusted R-squared=0.9489

其二：迴歸ANOVA表，雖然沒有作成表格形式，但上表的下方顯示實驗誤差(Residual standard error)=8.835，其平方就是MSe，誤差(迴歸稱殘差)自由度10，F值56.74，自由度迴歸項5、誤差項10，P-值為5.14e-07(非常顯著)。

其三：母迴歸係數的t檢定

另外殘差Residuals項其中位值為0.125，理論上應是0，所以有些偏差，對映Y值是不大的。

6 實驗結果的活用
ANOVA表與母迴歸係數顯著性檢定表已經將數據數理解析完成，確定A、C、D是關鍵因子，交互作用A×C與A×D顯著，一般先繪製主效果圖與交互作用圖(有顯著的)，方便決定最佳條件。
R有很多套件可以畫主效應圖與交互作用圖，本例為篩選設計所以我希望將4個因子畫在同一圖上，但又希望這些畫圖指令可運用在不同狀況，因此猶疑不決花了很多時間看套件與範例，最後還是沒有達到原本初衷，我的選擇是採用FrF2套件，其中MEPlot(主效果圖)與IAPlot(交互作用圖)二者適合於2水準的實驗資料畫圖，當然還有一些例如列出Alias Structure等重要功能。
如同第4節BsMD套件安裝，第一次執行MEPlot與IAPlot指令前需要事先安裝FrF2後執行指令

library(FrF2) MEPlot(lm.out) IAPlot(lm.out)

本來希望將MEPlot與IAPlot畫成一張圖但我沒做到，因此要注意執行MEPlot後應先將主效果圖複製或儲存然後再執行IAPlot，若二指令同時執行，主效果圖就被交互作用圖蓋住。
1) 主效果圖

此圖使用lm.out畫圖，lm.out是剔除不顯著項目B與一些交互作用的迴歸結果，因此只有A、C、D三個因子，當然也可以將四個因子全部畫出，做法是迴歸結果命名名稱不要一致(如lm.out與lm.out1)。

2) 交互作用圖

如同主效果圖，此交互作用圖也只有A、C、D三個因子的所有交互作用，因為A×C、A×D顯著所以最佳條件決定為A(1)C(-1)D(1)

7 最佳條件的推定

求最佳條件(1)C(-1)D(1)下的推定，可執行指令predict

其意為最佳條件(1)C(-1)D(1)的推定值為201.25，95%CI為(189.1956~213.3044)，此結果與Excel計算一致。

本例原為A~D四因子2水準完全因子設計共2⁴=16Runs，B因子不顯著而自模型剔除，因此成為2³×2，亦即三因子2水準重覆數2的完全因子實驗，每個實驗組合(Treatment)會有2個觀測值，本例A(1)C(-1)D(1)的2個觀測值是200與208，就習慣上
A(1)C(-1)D(1)的平均值為(200+208)/2=204，但母平均估計值是201.25，有別於觀測值的平均，此種母平均的推定值稱為Fit value(擬和值)或LS mean(最小平方的平均值)。

<全文玩>