如何用SAS軟件正確分析生物醫學科研資料Ⅺ.用SAS軟件實現正交設計定量資料的統計分析

柳偉偉，胡良平，陶麗新，毛瑋

1 正交設計的概念及應用場合

通常情況下，可認為正交設計所確定的因素間的水平組合是析因設計所確定的因素間的水平組合的一部分。正交設計就是用一系列規格化的正交表來安排各實驗因素及其水平組合的過程。與析因設計相比，正交設計是以不考察或少考察因素間高階交互作用項為代價換取很少實驗次數，從而達到安排多個實驗因素且取得較可靠實驗結果的目的。

在可以應用析因設計的實驗研究中，若所需要的實驗次數太多，并且高階交互作用可以忽略不計時，為了減少實驗次數，可以考慮選用正交設計。

2 正交表及其分類

正交表是正交設計的基本工具，它是根據均衡分散、整齊可比的思想經過嚴格的數學推導而編制出來的。正交表中的每一行代表各實驗因素的一種水平組合，稱為一個實驗點；正交表的每一列代表一種實驗效應，它可能代表某實驗因素、交互作用或實驗誤差的效應[1]。每張正交表都有一個表頭符號 Ln（Km），其中 L 代表正交表；n 代表正交表的行數，也就是實驗條件（或實驗點）數；m 代表正交表的列數，也就是正交表最多能安排的因素個數；K 代表正交表每一列中不同數字代碼的個數，也就是各因素的水平數。

正交表中總的自由度等于實驗次數減 1，如果每種實驗條件下只做一次獨立實驗，總的自由度就等于正交表的行數減 1；某一列的自由度等于該列的水平數減 1[2]。

一般來說，每張正交表都有一張與其對應的交互作用表，交互作用表用來說明任意兩列的交互作用列在正交表中所處的位置。

正交表可以分為同水平的正交表與混合水平的正交表兩大類。同水平正交表是指各因素的水平數相同，包括二水平正交表，例如 L4（23）、L8（27）、L16（215）；三水平正交表，例如 L9（34）、L27（313）；四水平正交表等等。混合水平正交表是指各因素的水平數不全相同，例如 L8（41×24）、L9（21×32）等等。實際應用時需要注意的是，有一類正交表沒有交互作用列，因此無法考察交互作用，例如 L12（211）、L18（37），它們屬于非標準的正交表。

3 正交設計的特點與實施步驟

正交設計有三個突出的特點，首先，由正交表挑出來的實驗點在空間上具有均勻分散性，也就是實驗點在實驗空間分布得很均勻，無論從哪個角度看，都是有代表性的點被挑選出來了；其次，由正交表挑出來的實驗點在統計分析時具有整齊可比性；第三，某些好的未包括在正交表中的實驗點，可以通過統計分析將其發現[1]。

在確定了實驗因素及其水平數之后，正交設計的實施可以分以下四個步驟來完成：第一步，根據因素數及各因素的水平數選擇合適的正交表，在能夠安排下所有實驗因素的基礎上，應該選擇最小號的正交表，使得需要進行的實驗次數較少。第二步，在選定的正交表的表頭上安排實驗因素及其交互作用，這一步稱為表頭設計，是所有步驟當中最關鍵的一步。如果因素間的交互作用可以忽略不計，各因素可以任意安排到各列中去。如果某些因素間存在交互作用，就需要根據交互作用表進行表頭設計。表頭設計的一個重要原則是避免混雜，混雜是指在正交表的同一列上安排了兩個或兩個以上的因素或交互作用[3]。當每種實驗條件下只做一次實驗時，如果要進行方差分析，就需要在正交表中安排至少一個空白列。在一些統計書上有現成的表頭設計表，可以根據這些表直接安排實驗。第三步，根據設計好的表頭，將標有單個實驗因素的那些列連同其下的“水平代碼”一起摘錄出來。第四步，結合實驗因素各水平的具體內容，將摘錄出來的各列的“水平代碼”轉換成“真實代碼”（即實驗因素的真實水平），并按正交表各行所決定的實驗條件進行具體實驗。

正交設計的實施過程，以及正交表與交互作用表將通過以下的實例進一步加以說明。

4 正交設計中一個值得注意的問題

人們常以為進行正交設計就不必要進行重復實驗了，即使正交表中未留下一個空白列時也是如此。這是過于迷信正交設計作用的表現。其實，通常（特指為節省樣本量時）正交設計就是析因設計的部分實施，換句話說，就是僅使用實驗因素全部水平組合當中的一部分，有時是很少的一部分，本質上屬于“對照不全”的實驗，正是由于正交設計方案挑選出來的那些實驗點具有“均勻分散、整齊可比”的特點，才使其具有較好的代表性，但并不意味著在結論的可靠性上它比全面實驗的析因設計還要好。析因設計中明確要求各實驗條件下至少要做 2 次獨立重復實驗，有兩個目的：第一，為了能夠比較準確地估計單個實驗因素的主效應和因素之間的各級交互作用的效應大小；第二，能夠比較真實地顯露各實驗條件下的實驗誤差的大小，提高結論的可信度。而正交設計實際上就是以犧牲第一點中的部分利益（即不估計高階交互作用，甚至部分低階交互作用也不能估計）為代價，換取較少的實驗次數，但它并沒有能力和資格承諾各實驗條件下僅做一次實驗，結果就一定是穩定的！也就是說，若在一個具體實驗中，僅當在特定實驗條件下所得到的實驗結果的波動是非常微小的（在專業上允許的范圍內，相當于確定性實驗結果，如在一個標準大氣壓下，水加熱到攝氏一百度時必然沸騰），才可不做重復實驗，否則必須做重復實驗。至于要做多少次重復實驗，取決于具體實驗的誤差大小和研究者對結果精確度的要求，最好根據一些基本數據和先驗知識，估計出重復實驗次數。

5 正交設計定量資料方差分析的SAS 實現

例1 在乙酰苯胺磺化工藝的研究中，有 4個實驗因素：反應溫度 A、反應時間 B、磺酸濃度 C、操作方法 D，各取兩水平。實驗的結果為產物的收率（%），希望弄清各因素在怎樣的搭配條件下收率最高。已知反應溫度與反應時間之間的交互作用不可忽視，各實驗條件下不必進行重復實驗，希望總實驗次數盡可能少一些。4個實驗因素的水平分別如下：

因素名稱（單位） 1 水平 2 水平反應溫度（℃） 50 70反應時間（h） 1 2磺酸濃度（%） 17 27操作方法 攪拌 不攪拌

正交設計的實施：本例如果選擇析因設計，不同的實驗條件數為24=16 種，各實驗條件下至少要做 2 次獨立重復實驗，總實驗次數至少為32 次，實驗次數相對較多，與要求不符。由于只需要考慮兩個實驗因素之間的一階交互作用，高階交互作用可以忽略，故選擇正交設計比較合適。

研究中的4個實驗因素都是兩水平，應該選擇二水平正交表。由于只需要考慮兩個實驗因素之間的一階交互作用，連同 4個兩水平因素，共需要安排 5 項，每項的自由度都是 1，將占用 5個自由度。在二水平正交表中，每列的自由度為1，故需要 5 列來安排實驗因素及其交互作用。如果選擇 L4（23）正交表不妥，因為該正交表只有 3 列，自由度為3，用來安排該實驗是不夠的；若選擇 L16（215）正交表也不妥，因為該正交表有 15 列，自由度為15，遠遠大于所需要的5 列，意味著實驗次數將較多；若選擇 L8（27）正交表，是非常合適的，因為該正交表有 7 列，自由度為7，比所需要的5 列多出兩列，正交表中將有兩個空列，可用于估計實驗誤差，此時，只需要做8次實驗，就可以滿足研究者的要求。L8（27）正交表是符合要求的最小號的正交表，具體見表1，其交互作用表見表2。

表1 L8（27）正交表

表2 L8（27）正交表的交互作用表

L8（27）正交表有 8 行 7 列，用該表來安排實驗，要在8 種實驗條件下總共至少做 8 次實驗，最多可以安排的因素個數為7，每個因素都有兩個水平，分別用 1、2 來表示。

在表 2 中，第 1個列號是由上到下編號的，代表著 7行；第 2個列號是從左至右編號的，代表 7 列。無論是 7行還是 7 列指的都是 L8（27）正交表中的“7 列”，也就是說，把 L8（27）表中的7個列號同時放在橫向與縱向兩個方向上。如果要查任意兩列的交互作用列，只需要找到這兩列列號的交叉位置的數字就可以了。例如，查第 6 列和第 7 列的交互作用列，可以橫向看第 6 行、縱向看第 7列，其交叉位置上的數為1，說明這兩列的交互作用列為第1 列。值得注意的是，交互作用項的自由度為兩個因素的自由度之積，每個二水平因素的自由度都為1，因此兩個二水平因素的交互作用項僅占 1 列；兩個三水平因素的交互作用項的自由度為2×2=4，因此兩個三水平因素的交互作用項需占 2 列；同理，兩個四水平因素的交互作用項需占用 3 列，依此類推[1]。

根據表1和表 2，就可以進行表頭設計了。將反應溫度 A 放在L8（27）表的第 1 列，反應時間 B 放在第 2 列。然后查交互作用表，可知反應溫度與反應時間的交互作用AB 應落在第 3 列，因此第 3 列不能再安排其他因素，否則會產生混雜。于是可把磺酸濃度 C 放在第 4 列，操作方法 D 可以放在5、6、7 三列中任何一列上，不妨將 D 放在第 7 列上。這樣就完成了該實驗的表頭設計，如表 3所示。

表3 用 L8（27）正交表安排乙酰苯胺磺化實驗的表頭設計

完成表頭設計之后，將 L8（27）表中的第 1、2、4、7列摘錄出來，將表中的代碼 1、2 轉化成 4個實驗因素的真實水平，然后按照表中的實驗點安排實驗，得到實驗結果見表4。

表4 L8（27）正交設計的實驗結果

對該資料進行分析的SAS 程序如下：

data prg1;input a b c d y@@;cards;111165 112274121271 122173211270 212173221162 222267;run;proc glm data=prg1;class a b c d;model y=a b c d a*b;run;

程序說明：數據步建立數據集 prg1，變量 a、b、c、d、y 分別代表反應溫度、反應時間、磺酸濃度、操作方法和結果變量收率。這里各個因素的取值是用正交表中的水平代碼1、2 表示的，這樣表示的優點是寫起來比較簡單。當然也可以采用各個因素的真實水平，只不過輸入時會繁瑣一些。需要注意的是，交互作用列和空白列中的數字是不出現在數據集中的。與其他各種設計類型的定量資料一樣，進行正交設計定量資料的方差分析時仍然采用 GLM 過程。在class語句中，指定分組變量為a、b、c、d。使用 model 語句指定模型的具體形式，該語句等號左端為結果變量 y，等號右端為需要分析的效應，包括四個主效應和一個交互效應。

主要輸出結果與結果解釋：

The GLM procedure dependent variable: y

以上是對整個模型進行假設檢驗的結果，輸出內容包括自由度、離均差平方和、均方、F 值和P 值，其中 F=11.82、P=0.0798。

R-square Coeff var Root MSE y Mean 0.967276 2.101244 1.457738 69.37500

以上是一些描述性統計量，R-square 為決定系數，Coeff var 為變異系數，Root MSE 為誤差均方的平方根，y Mean為結果變量的均數。

Source DF Type III SS Mean square F value Pr ＞F a 1 15.12500000 15.12500000 7.12 0.1165 b 1 10.12500000 10.12500000 4.76 0.1607 c 1 45.12500000 45.12500000 21.24 0.0440 d 1 10.12500000 10.12500000 4.76 0.1607 a*b 1 45.12500000 45.12500000 21.24 0.0440

以上是對各個效應進行假設檢驗的結果，其中因素 A對應的F=7.12、P=0.1165；因素 B和因素 D 對應的F=4.76、P=0.1607。這三者對于結果變量的作用沒有統計學意義。因素 C 與交互作用 AB 對應的F=21.24、P=0.0440＜0.05，這兩項對結果變量的作用有統計學意義。根據上述結果可知，磺酸濃度、反應溫度與反應時間的交互作用對產物的收率存在影響，而反應溫度、反應時間及操作方法對收率沒有影響。還可對 P ＞0.05的項進行逐一淘汰，通常從最大 P 值對應的項開始淘汰，直到所有項對應的P 值都小于0.05 為止。

對正交設計的定量資料除了進行方差分析外，還可以通過簡單的計算進行直觀分析，直觀分析包括綜合比較和極差分析。綜合比較是指比較不同因素各水平的結果大小，篩選全部因素水平最佳組合條件，此時可不必做復雜的方差分析。一個因素的極差是指該因素各水平均值的最大值與最小值之差，極差值越大，說明改變這一因素的水平會使結果產生較大的變化[2]。極差分析就是比較各因素極差的大小，可以區分因素的主次。現將例1 中不同因素的兩個水平的實驗結果分別求和，具體見表5。

表5 L8（27）正交設計的直觀分析計算表

表5 中 T1m為第 m 列水平數為1 時實驗結果的合計，T2m為第 m 列水平數為2 時實驗結果的合計。例如，反應溫度 A 為1 水平，也就是 50 ℃ 時實驗結果的合計為283，其為2 水平時實驗結果的合計為272。

根據表 5 中的數據進行綜合比較，由于T11（283）＞T21（272），說明反應溫度 50 ℃ 比 70 ℃ 時的產物收率高；T12（282）＞T22（273），說明反應時間 1 h 比 2 h的產物收率高；T13（268）＜T23（287），說明反應溫度和反應時間存在交互作用；T14（268）＜T24（287），說明磺酸濃度 27%比 17% 時的產物收率高；T17（273）＜T27（282），說明不攪拌比攪拌時的產物收率高。單因素分析的結論為：反應溫度 50 ℃、反應時間 1 h、磺酸濃度 27%、操作方式為不攪拌時收率較高。但是因為反應溫度和反應時間存在交互作用，還需再按表 5 第 1、2 兩列計算 A、B 兩個因素不同水平組合下收率的合計：

反應溫度A反應時間 B 50 ℃ 70 ℃1 h 65 + 74=139 70 + 73=1432 h 71 + 73=144 62 + 67=129

4 種水平組合下的合計結果顯示，反應溫度 50 ℃ 與反應時間 2 h的收率最高，這與單因素分析的結論并不一致，所以反應溫度 50 ℃、反應時間 2 h、磺酸濃度 27%、操作方式為不攪拌是最佳組合條件[4]。同時也可以看出，這個組合條件在本次研究中并沒有安排進行實驗。

由表 5 可以算得 A、B、C、D 這 4個實驗因素的極差分別為RA=2.75、RB=2.25、RC=4.75、RD=2.25，RC＞RA＞RB=RD，所以磺酸濃度是主要因素，反應溫度次之，反應時間和攪拌方式是最次要的因素。

例2 在一項微生物培養液成分的優化實驗中，希望找出最優的培養液成份，欲考慮的因素及其水平見表6。已知交互作用 AC 存在的可能性極大，AB和AE 存在的可能性不大，但無把握斷定其不存在，希望通過實驗加以考察[5]。

表6 培養液優化實驗的因素水平表

表7 L27（313）正交表

正交設計的實施：本實驗涉及5個三水平的實驗因素，若進行析因設計，則至少需要做 35×2=486 次實驗，實驗次數過多。由于只需要考慮 3個一階交互作用，故適合采用正交設計。

本例中每個實驗因素的自由度都為2，5個因素的自由度為10；每個交互作用的自由度為2×2=4，3個一階交互作用的自由度為12；總的自由度為22。在三水平正交表中，每一列的自由度為2，故需要 11 列來安排實驗因素及其交互作用。采用 L27（313）正交表安排此實驗是可行的，因為此表共有 13 列，5個因素占用 5 列，3個一階交互作用占用 6 列，還剩兩個空列用于估計實驗誤差，具體見表7。同時，假設已經得到了現有的L27（313）正交表的表頭設計表，具體見表8。

表8 L27（313）正交表的表頭設計表

在表 8 中，每個實驗因素占 1 列，每個一階交互作用占 2 列，例如因素 A 與B的交互作用 AB 占用第 3 列和第 4 列，分別用 AB1和AB2表示。根據此表頭設計，可以直接將各因素安排到相應的列中，然后摘錄出單個實驗因素對應的列，將各列的水平代碼轉換成實驗因素的真實水平。限于篇幅，經轉換后的表格不再給出，僅將實驗結果列于表 7 中。

對該資料進行分析的SAS 程序如下：

data prg2;input a b c d e y@@;cards;111110.69112220.54113330.37121320.66122130.75123210.48131230.81132310.68133120.39211110.93212221.15213330.90221320.86222130.97223211.17231230.99232311.13233120.80311110.69312221.10313330.91321320.86322131.16323211.30331230.66332311.38333120.73;run;proc glm data=prg2;class a b c d e;model y=a b c d e a*c a*b a*e;run;proc glm data=prg2;class a b c d e;model y=a b c d e a*c;run;

程序說明：數據步建立數據集 prg2，與例1 一樣，這里各個因素的取值仍然采用正交表中的水平代碼 1、2、3 表示。在L27（313）正交表中，因素 E和因素 D 分別位于第 8 列與第 11 列，E的位置在D 之前，所以輸入數據時要注意它們的順序。本例使用了兩個 GLM 過程步，兩者的不同之處在于第二個 GLM 過程去掉了交互作用項AB 與AE。

主要輸出結果與結果解釋：

Source DF Sum of squares Mean square F value Pr ＞F Model 22 1.75077037 0.07958047 4.97 0.0651 Error 4 0.06408148 0.01602037 Corrected total 26 1.81485185 Source DF Type III SS Mean square F value Pr ＞F a 2 0.89516296 0.44758148 27.94 0.0045 b 2 0.05031852 0.02515926 1.57 0.3138 c 2 0.23000741 0.11500370 7.18 0.0475 d 2 0.06667407 0.03333704 2.08 0.2402 e 2 0.10738519 0.05369259 3.35 0.1397 a*c 4 0.30961481 0.07740370 4.83 0.0781 a*b 4 0.04577037 0.01144259 0.71 0.6239 a*e 4 0.04583704 0.01145926 0.72 0.6233

以上是第一個 GLM 過程的主要輸出結果，對整個模型進行假設檢驗的F=4.97、P=0.0651。在各個效應項的檢驗結果中，因素 A 對應的F=27.94、P=0.0045；因素 C對應的F=7.18、P=0.0475。這兩者的作用都有統計學意義。交互作用 AC的P 值接近于0.05的臨界水平。交互作用 AB 與AE 沒有統計學意義，結合實驗之前已經掌握的情況，將 AB 與AE 去掉之后重新擬合模型。

Source DF Sum of squares Mean square F value Pr ＞F Model 14 1.65916296 0.11851164 9.13 0.0002 Error 12 0.15568889 0.01297407 Corrected total 26 1.81485185 Source DF Type III SS Mean Square F value Pr ＞F a 2 0.89516296 0.44758148 34.50 ＜0.0001 b 2 0.05031852 0.02515926 1.94 0.1863 c 2 0.23000741 0.11500370 8.86 0.0043 d 2 0.06667407 0.03333704 2.57 0.1178 e 2 0.10738519 0.05369259 4.14 0.0430 a*c 4 0.30961481 0.07740370 5.97 0.0070

以上是第二個 GLM 過程的輸出結果，對整個模型進行假設檢驗的F=9.13、P=0.0002。在對各個效應的檢驗中，因素 A、C、E 以及交互作用 AC 都有統計學意義，其 P 值分別為＜0.0001、0.0043、0.0430、0.0070；因素 B與D的作用沒有統計學意義。最好，將 B、D 淘汰后再進行分析，以得到更穩定的計算結果，此處從略。

[1]Hu LP.Application of statistical triple-type theory in the experiment design.Beijing: People’s Military Medical Press, 2006:121-132.(in Chinese)胡良平.統計學三型理論在實驗設計中的應用.北京:人民軍醫出版社, 2006:121-132.

[2]Mao SS, Zhou JX, Chen Y.Design of experiment.Beijing: China Statistics Press, 2004:117-200.(in Chinese)茆詩松,周紀薌, 陳穎.試驗設計.北京: 中國統計出版社, 2004:117-200.

[3]Ren LQ.Optimum design and analysis of experiments.2nd ed.Beijing: Higher Education Press, 2003:10-58.(in Chinese)任露泉.試驗優化設計與分析.2版.北京: 高等教育出版社, 2003:10-58.

[4]Sun ZQ.Medical statistics.Beijing: People’s Medical Publishing House, 2002:188-194.(in Chinese)孫振球.醫學統計學.北京:人民衛生出版社, 2002:188-194.

[5]Wang WZ, Mao SS, Zeng LR.Design and analysis of experiments.Beijing: Higher Education Press, 2004:114-191.(in Chinese)王萬中, 茆詩松, 曾林蕊.試驗的設計與分析.北京: 高等教育出版社, 2004:114-191.