線性回歸分析法在正交試驗中的應(yīng)用

范步高

(上海長征富民金山制藥有限公司，上海201506)

正交試驗數(shù)據(jù)處理的常用方法一般有極差分析 (直觀分析)和方差分析。極差分析簡便、直觀，但不能區(qū)分指標的差異究竟是因子的水平變化所引起還是由試驗的誤差所引起，故極差分析的準確性較低，所以，為提高正交試驗結(jié)果分析的準確性應(yīng)首選方差分析。但正交試驗的方差分析一般須預(yù)留空白列 (誤差列和交互作用列)，代價是試驗總數(shù)增加或減少獨立因子數(shù)，否則，正交試驗方差分析的準確性也將大大降低，此時較為理想的選擇則應(yīng)是“線性回歸分析法”［1-4］。

1 實例分析

《中成藥》2011年 (第33卷)第7期第1206頁載文“正交設(shè)計法優(yōu)選苦參炮制工藝的研究”［5］(以下簡稱原文)，原文目的是優(yōu)選苦參的切制工藝，選用L9(34)正交表進行優(yōu)選，現(xiàn)將原文的表2(因素和水平考察)、表3(苦參炮制工藝正交試驗結(jié)果)和表4(苦參炮制工藝方差分析)合并為本文的表1，原文直觀分析 (極差分析)與方差分析得出的最佳條件與直接觀察法相同，均為A1B3C3，即第3號試驗，Y3=Ymax=1.831，因素重要性排序為A＞C＞B。3批驗證結(jié)果:苦參堿和氧化苦參堿指標分別為 1.83%、1.80%、1.82%，均值為 1.82%，RSD為0.64%。

1.1 方差分析的正確計算 原文為3因子3水平，選擇L9(34)正交表，前提條件應(yīng)是因子間無兩兩交互作用或兩兩交互作用可以忽略不計，另外3因子排畢還應(yīng)剩余“空列”1列，可以作為誤差項進行方差分析 (原文將“空列”遺漏，故計算有誤)，因B和C因子的均方差均小于“空列”，為防止誤判因子的顯著性，提高方差分析的準確性，故合并三列為誤差項 (SS空+SSB+SSC=SSe)。結(jié)果見表2。

表1 原文苦參炮制工藝正交試驗結(jié)果分析

表2 苦參炮制工藝正交試驗方差分析

由表2可知，因素A、B、C水平的變化對指標的影響均不顯著 (與原文結(jié)論差異較大)，即差異較多是由試驗誤差過大所引起。所以，對不顯著因子的水平優(yōu)選一般不是選擇相對較高指標值所對應(yīng)的水平，而是綜合考慮該水平在實際生產(chǎn)中的生產(chǎn)周期、能耗、質(zhì)量可控性等情況而選定。所以，若綜合考慮苦參炮制工藝的最佳條件應(yīng)首選A1B1C1(浸泡20 min，軟潤16 h，60℃干燥)，即該輪L9(34)正交試驗中第1號試驗，Y1=1.732，低于第3號試驗指標值Y3=1.831，但此結(jié)果僅僅是在因子間無兩兩交互作用或兩兩交互作用可以忽略的前提下。根據(jù)L9(34)正交表任意兩列間的交互作用列為另外兩列，現(xiàn)存在MS空＞MSC＞MSB，故原文前提條件不成立，即 A、B、C因子的兩兩交互作用不可忽略。

表2結(jié)果可由表1原始數(shù)據(jù)用Excel統(tǒng)計函數(shù)計算得到［6］，也可仿照原文由簡單公式計算:T=ΣYi=26.512，CT=T2/n=26.447;總平方和SST=T2－CT;A、B、C及空列(第4列)平方和SSjCT。對規(guī)格化正交表存在等式:SST=Σ SSj，可以用于檢驗計算的正確與否。

1.2 線性回歸分析及預(yù)測 對三水平的A、B、C因子間的交互作用，解決的辦法可以重新選擇較大的L18(37)或L27(313)正交表再進行一次有交互作用的正交試驗進行優(yōu)選，也可以在原L9(34)正交試驗的基礎(chǔ)上，利用多元線性回歸分析法進行工藝優(yōu)選和結(jié)果預(yù)測。

1.2.1 多元線性回歸分析模型構(gòu)建 根據(jù)正交試驗的方差分析結(jié)果，假設(shè)A、B、C存在較強的兩兩交互作用AB、AC和 BC，回歸子集大小為六元 (A、B、C、AB、AC、BC)，線性回歸模型為Y*=b0+b1A+b2B+b3C+b4AB+b5AC+b6BC，經(jīng)用矩陣法求得各因子的偏決定系數(shù) (偏相關(guān)系數(shù)的平方)進行顯著性檢驗和篩選，“最優(yōu)”子集為五元 (A、B、C、AB、AC)、“次優(yōu)”子集為四元 (B、C、AB、AC)(表3，表4)。

表3 線性回歸“最優(yōu)”子集Excel篩選模型

(續(xù)表3)

表4 線性回歸“最優(yōu)”子集篩選結(jié)果

1.2.2 線性回歸分析結(jié)果及置信區(qū)間 由表3可得五元“最優(yōu)”線性回歸方程為 Y*=1.209 6+0.018 6A+0.095 2B－0.010 0C－0.003 0AB+0.000 3AC，回歸標準誤差Sy=0.014 5，回歸方程顯著性P值＜0.01，各回歸因子顯著性Pj值＜0.05，因子重要性排序AB＞B＞C＞AC＞A。對原文優(yōu)選水平A1B3C3進行預(yù)測，均值落在其預(yù)測區(qū)間內(nèi)(表3)，證明回歸方程的預(yù)報可信。參考方差分析結(jié)果，經(jīng)多次單因素試驗回歸得優(yōu)選水平為A1B3C1(浸泡20 min，軟潤20 h，60℃干燥)，此水平組合未包含在L9(34)正交試驗中，故需進行至少3批的驗證試驗加以確定，其驗證試驗結(jié)果預(yù)測值約為2.026%，95%置信區(qū)間在1.997%～2.054%，高于第3號正交試驗 (A1B3C3)最高值Y3=1.831%(見表3)。

1.2.3 表3中Excel函數(shù)或公式設(shè)置

(1)試驗因子及交互項水平數(shù)據(jù)與指標結(jié)果:A、B、C及Y數(shù)據(jù)可直接輸入，交互作用的因子數(shù)據(jù)利用乘法公式輸入，E3=B3×C3(意即E3單元格輸入“=B3×C3”，下同)，F(xiàn)3=B3×D3，G3=C3×D3，分別確認后選中E3∶G3，向↓拖曳復(fù)制至G11。

(2)相關(guān)系數(shù)矩陣1:B14=CORREL(S|BS|3∶S|BS|11，BS|3∶BS|11)，確認后選中并向↓拖曳復(fù)制至B20，并從B15起依次將S|BS|3∶S|BS|11中的列號B分別修改為C、D、E、F、G、H，分別確認后選中B14∶B20，向→拖曳復(fù)制至H列。

(3)相關(guān)系數(shù)逆矩陣1:選中B22∶H28(區(qū)域大小須與B14∶H20矩陣相同)，輸入:=MINVERSE(B14:H20)，三鍵確認 (一般特定公式或數(shù)組公式輸入均需shift+ctrl+enter三鍵確認，當判別有疑時也可用三鍵確認解決，下同)。

(4)回歸因子顯著性檢驗1:偏決定系數(shù)B30=B28^2/B22/S|H28(若求偏相關(guān)系數(shù)，B30=－B28/SQRT(B22/S|H28))，偏回歸因子F值B31=B30/(1－B30) × (9－S|H31－1)，偏回歸因子 P值 B32=FDIST(B31，1，9－S|H31－1)，偏回歸因子重要性排序 B33=RANK(B30，S|B30:S|G30)，分別確認后選中B30∶B33，向→拖曳復(fù)制至G列。Mpi統(tǒng)計量 H33=SUM(B32:G32)/H31× (9+H31+1)。

(5)相關(guān)系數(shù)矩陣2:根據(jù)回歸因子顯著性檢驗1的結(jié)果，去掉相關(guān)系數(shù)矩陣1中最不顯著因子BC(第六行及第六列數(shù)據(jù))，利用等號將剩余標志和數(shù)據(jù)“移植”，如C35=B13，C41=B20。為方便回歸因子顯著性檢驗公式無需修改地重復(fù)使用，需將原相關(guān)系數(shù)矩陣1中第七行數(shù)據(jù)上移至第六行、原第1～5列數(shù)據(jù)整體右移。

(6)相關(guān)系數(shù)逆矩陣2:選中C43∶H48，輸入=MINVERSE(C36:H41)，三鍵確認。

(7)回歸因子顯著性檢驗2:復(fù)制表3(4)，清除B50∶B53，修改子集m大小，H53復(fù)制H33后粘貼。

(8)五元“最優(yōu)”子集回歸統(tǒng)計量:選中5行 (m+1)列即 C56∶H60，輸入:=LINEST(S|HS|3∶S|HS|11，S|BS|3∶S|FS|11，TRUE，TRUE)，三鍵確認?；貧w方程P值:C61=FDIST(C59，9-1-D59，D59)

(9)試驗優(yōu)水平篩選及預(yù)測:利用多次單因素試驗進行優(yōu)水平篩選，并假設(shè)殘差服從正態(tài)分布，U0.05=1.96，則95%置信度預(yù)測值誤差=1.96×Sy，所以，預(yù)測值F64=SUM(S|CS|56∶S|GS|56×A64∶E64)+S|HS|56，三鍵確認;95%下限值G64=S|F64－1.96×S|DS|58，95%上限值H64=S|F64+1.96×S|DS|58，分別確認后選中F64∶H64，并向↓拖曳復(fù)制至行72。

對于確定的“最優(yōu)”子集，偏回歸系數(shù)等回歸統(tǒng)計值也可以直接利用Excel回歸分析工具自動求出，可以作為對以上顯著性檢驗計算正確與否的檢驗 (自動給出的偏回歸系數(shù)t檢驗與回歸因子的偏決定系數(shù)的F檢驗等值)，方法是:分別點擊Excel工具→數(shù)據(jù)分析→回歸分析→確定，填寫回歸對話框，Y值區(qū)域輸入S|BS|1∶S|BS|10，X值區(qū)域輸入S|CS|1∶S|GS|10，勾選“標志”，點選“輸出區(qū)域”，點擊B73，點擊確定，即可得到更為詳細的回歸分析結(jié)果 (略)。

2 討論

2.1 工藝優(yōu)化水平的確定 對q個水平m*個獨立回歸因子的線性回歸方程進行優(yōu)化工藝水平的確定，既可采用qm*次單因素試驗法，也可采用更為準確的qm*次全面試驗法或“規(guī)劃求解法”［3］，也可以采用“極值法”{4}或根據(jù)方程的復(fù)雜程度靈活組合應(yīng)用，如五元“最優(yōu)”線性回歸方程Y*=1.209 6+0.018 6A+0.095 2B－0.010 0C－0.003 0 AB+0.000 3AC，若求其極值較為不易，但在試驗水平考察范圍內(nèi)有A↘Y*↗ (Y*對A的一階偏導(dǎo)小于零)、C↘Y*↗ (Y*對C的一階偏導(dǎo)小于零)和B↗Y*↗ (Y*對B的一階偏導(dǎo)大于零)的單調(diào)性趨勢［9］，所以，Y*最大值所對應(yīng)的優(yōu)化水平可選定為A1B3C1，這也與一般中藥的穩(wěn)定性規(guī)律 (低溫干燥有利于減少氧化、分解等反應(yīng)對中藥有效成分群所造成的損失)和苦參堿及氧化苦參堿的理化性質(zhì) (極性大，易水溶，易水解)相吻合。

多次單因素試驗法是先固定 (m*－1)個獨立因子的水平并與另一因子的各個水平進行q次全搭配預(yù)測，根據(jù)預(yù)測結(jié)果得出該因子的較優(yōu)水平并固定，然后更換下一個獨立因子重復(fù)上述全搭配預(yù)測，直至進行qm*次 (包括m*－1次重復(fù))得出所有m*獨立因子的較優(yōu)水平搭配即為回歸分析法的優(yōu)選水平組合，但當回歸方程在水平考察范圍內(nèi)存在拐點［9］時，不同水平的固定將會影響優(yōu)選水平的準確性。

2.2 “最優(yōu)”線性回歸方程的篩選準則 “最優(yōu)”線性回歸方程的篩選準則通常有復(fù)決定系數(shù)R2、復(fù)相關(guān)系數(shù)R、回歸F值等愈大愈好，子集m、殘差平方和SSe、平均殘差平方和 MESS(均方 MSe)、平均預(yù)測均分誤差 MPESS、AIC準則、BIC準則等愈小愈好［7-8］。本文的篩選方法或步驟:①選擇“最優(yōu)”子集m，使m子集內(nèi)的各個因子的顯著性Pj值均小于0.05或0.10，②根據(jù)m“少而精”之原則，設(shè)定統(tǒng)計量Mpi=(n+m+1)×ΣPj/m，選擇Mpi最小子集;③選擇線性回歸方程的顯著性P值小于0.05的“最優(yōu)”子集 (各因子顯著性Pj值均小于0.05的“最優(yōu)”子集，其對應(yīng)的線性回歸方程的顯著性P值一般也小于0.05，反之，則可能性降低，需要進一步對各因子進行顯著性檢驗加以確定)。

對于三獨立因子 (m*=3)試驗，其二次多項式回歸因子 (一次項、二次項、交互項)共有 9項，若一次項A、B、C和交互項AB、AC、BC子集中未能在表3模型中篩選出“最優(yōu)”子集，可以將二次項A2、B2、C2回歸因子替換表3相關(guān)系數(shù)矩陣1中最小偏決定系數(shù)對應(yīng)的回歸因子或其他因子 (可以重復(fù)使用復(fù)制、粘貼、撤銷鍵，表中公式無需修改)進行篩選。

回歸因子m要求“少而精”，實則是要求“最優(yōu)”線性回歸方程內(nèi)的自變量因子均為顯著因子，而方程外的所有因子均為不顯著因子，同時也可以認為是要求線性回歸方程的殘差自由度 (離回歸自由度)盡可能的大 (dfe=n－m－1)，目的是避免擬合的方程回歸不錯 (殘差較小)，預(yù)報不好 (驗證試驗結(jié)果誤差較大)，其解決的辦法是增加觀察值 (試驗次數(shù))以達到n≥2m或至少保證殘差自由度(n－m－1)≥2。

2.3 正交試驗數(shù)據(jù)處理方法的選擇

正交試驗的目的，①是考察各因子的水平差異對指標值影響的相對大小，②是篩選最優(yōu)化試驗條件。

極差分析計算最為簡單，可達①之目的，但極差分析未及消除試驗誤差對指標值的影響，故其極差較大的主要因子未必是顯著因子［10］。所以，欲達②和①之目的，經(jīng)典的方差分析則較為準確，應(yīng)為首選。

在進行正交試驗的Ln(qm)正交表中，獨立因子各占1列，交互因子各占 (q－1)d列 (一般僅考慮兩兩交互作用，其交互作用級數(shù)d=2－1)，誤差項也至少須占1列(空列)，根據(jù)n(試驗總次數(shù)或正交表行數(shù))、q(因子水平數(shù))和m(因子列數(shù)或正交表列數(shù))三者存在n=qk或n=sq2(k=2，3，4，…，s=1，2，3，…)或n=m(q－1)+1的關(guān)系式［11］可知，m和q增加，n將呈番數(shù)或指數(shù)增長，結(jié)果是造成試驗成本大大增加。所以，正交試驗的方差分析也僅在因子間沒有交互作用或交互作用可以忽略，或者獨立因子數(shù)和水平數(shù)均較少時較為實用。

線性回歸分析與方差分析都是研究數(shù)據(jù)的統(tǒng)計方法，線性回歸法常用于均勻設(shè)計，方差分析法則常用于正交設(shè)計，但當正交設(shè)計的獨立因子 (自變量，非隨機變量)為連續(xù)性變量且與指標變量 (因變量，隨機變量)存在顯著線性關(guān)系時也可以采用線性回歸分析法［12］，其優(yōu)點之一是正交表的m列可以排滿獨立因子而達包括交互因子優(yōu)選水平在內(nèi)的②和①之目的，二是根據(jù)正交試驗信息量可以較為方便的構(gòu)建線性回歸模型，用于工藝過程預(yù)測和控制，并能指導(dǎo)超越正交試驗原有水平的考察范圍捕捉到更優(yōu)化的工藝條件，如上述優(yōu)選水平A1B3C1，均為“邊界”水平組合，在進行驗證試驗時，可以根據(jù)實踐經(jīng)驗同時考慮選擇小于A1(如A0:浸泡10min)或C1(如C0:55℃真空干燥)或大于B3(如B4:軟潤24 h)的水平組合進行對照。經(jīng)預(yù)測，A0B4C0苦參炮制工藝更優(yōu)，回歸預(yù)測值為2.570%，95%置信區(qū)間在2.542% ～2.599%，高于原優(yōu)選水平A1B3C1的預(yù)測值2.026%(表3)，這也是線性回歸分析優(yōu)于正交試驗方差分析的一個結(jié)果體現(xiàn)。

［1］李云雁，胡傳榮．試驗設(shè)計與數(shù)據(jù)處理［M］．北京:化學(xué)工業(yè)出版社，2005．

［2］胡譽滿，謝曉鳴．利用回歸分析對正交試驗試驗結(jié)果進行修正［J］．工科數(shù)學(xué)，2000，16(4):32-34．

［3］俞鐘行．質(zhì)量工程師教材例2.3-1新解［J］．質(zhì)量春秋，2011(11):48-351．

［4］郭俊旺，劉曉峰，魏彩云．正交設(shè)計L9(34)的漸進優(yōu)化線性回歸方法［J］．計算機與應(yīng)用化學(xué)，2010，27(11):1503-1508．

［5］鄧捷圓，胡 馨，張英華，等．正交設(shè)計法優(yōu)選苦參炮制工藝的研究［J］．中成藥，2011，33(7):1206-1208．

［6］范步高．正交試驗方差分析的Excel通用計算與應(yīng)用［J］．中國醫(yī)藥工業(yè)雜志，2011，42(10):793-795．

［7］李東風(fēng)，鄭忠國．最優(yōu)線性回歸的計算方法［J］．數(shù)理統(tǒng)計與管理，2008，27(1):87-95．

［8］丘冠英．“最優(yōu)”回歸方程的選擇準則和具體方法［J］．宜春學(xué)院學(xué)報:自然科學(xué)，2003，25(6):26-27．

［9］王順鳳，夏大峰，朱鳳琴，等．高等數(shù)學(xué)(上)［M］．北京:清華大學(xué)出版社，2009．

［10］郝拉娣，張 嫻，劉 琳．科技論文中正交試驗結(jié)果分析方法的使用［J］．編輯學(xué)報，2007，19(5):340-341．

［11］全國質(zhì)量專業(yè)技術(shù)人員職業(yè)資格考試辦公室．質(zhì)量專業(yè)理論與實務(wù)(中級)［M］．北京:中國人事出版社，2006．

［12］喬克林，呂 佳．方差分析與回歸分析之比較［J］．延安大學(xué)學(xué)報:自然科學(xué)版，2009，28(2):34-36．