0802 回歸分析簡介及其在均勻設(shè)計中的應(yīng)用

1、第二節(jié)回歸分析簡介及其在均勻設(shè)計中的應(yīng)用回歸分析是數(shù)據(jù)分析的有力工具，它能揭示變量之間的相互關(guān)系，因此在均勻設(shè)計的數(shù)據(jù)分析中成為主要的手段，回歸分析方法和理論十分豐富，有關(guān)書籍數(shù)以百計，這里僅作一梗概介紹，細節(jié)可以參看有關(guān)書籍，如26，29，30數(shù)據(jù)處理可使用統(tǒng)計軟件包SAS，SPSS，MINITAB，BMDP，S等，國內(nèi)許多部門如中國均勻設(shè)計學(xué)會為均勻設(shè)計及其數(shù)據(jù)分析制作了專用統(tǒng)計軟件包，使用更為方便。2.1 一元線性回歸模型 由于均勻設(shè)計的數(shù)據(jù)分析要利用回歸分析，因此需要對回歸分析作一扼要介紹。一元線性回歸是處理兩個變量之間關(guān)系的最簡單的模型。本章將詳細討論這個模型。一元線性回歸雖簡單，但

2、從中可以了解回歸分析方法的基本思想/方法和應(yīng)用。 我們首先通過一個例子說明如何建立一元線性回歸方程。 例3 為了估計山上積雪融化后對下游灌溉的影響，在山上建立了一個觀測站，測量了最大積雪深度(X)與當年灌溉面積(Y)，得到連續(xù)10年的數(shù)據(jù)。 年序最大積雪深度X(尺)灌溉面積Y(千畝)115.228.6210.419.3321.240.5418.635.6526.448.9623.445.0713.529.2816.734.1924.046.71019.137.4 為了研究這些數(shù)據(jù)中所蘊含的規(guī)律性，我們把各年最大積雪深度作橫坐標，相應(yīng)的灌溉面積作縱坐標，將這些數(shù)據(jù)點標在平面直角坐標圖上，如圖9，

3、這個圖稱為散點圖。 從圖9看到，數(shù)據(jù)點大致落在一條直線附近，這告訴我們變量X與Y之間的關(guān)系大致可看作是線性關(guān)系，從圖9還看到，這些點又不都在一條直線上，這表明X與Y的關(guān)系并沒有確切到給定X就可以唯一地確定Y的程度。事實上，還有許多其他因素對Y產(chǎn)生影響，如當年的平均氣溫，當年的降雨量等等，這些都是影響Y取什么值的隨機因素。如果我們只研究X與Y的關(guān)系，可以假定 有如下結(jié)構(gòu)式：Y=+X+ (2.1) 式中， 稱為回歸系數(shù)，X為自變量，Y為因變量，表示隨機誤差，常常假定遵從正態(tài)分布N(0，2)，這表示誤差為正和負的機會一樣多，2 表示誤差的大小。式中，2 通常是未知的，它們要通過數(shù)據(jù)的信息來估計。 設(shè)

4、(Xi,Yi)，i=1，n為一組數(shù)據(jù)，若用回歸方程(2.1)來擬合，則當X=時的估計值為 (2.2)自然，我們希望求和使與很接近.也就是說，我們要決定一條直線，使其與所有的點都比較接近，最流行求， 估計值的辦法是用最小二乘法，令 (2.3)最小二乘法是求和使Q達極小，使Q達極小的和值記為a和b.利用微積分中求極值的辦法求得 (2.4)式中 (2.5)利用這些公式到例3，得:于是 b=415.606/230.656=1.802×18.88=2.511從而回歸方程為 試將該直線畫在圖9上，可以看到擬合的效果是不錯的，衡量擬合效果的好壞，如下的方法是十分有用的。 (a) 相關(guān)系數(shù) 相關(guān)系數(shù)

5、 用于描敘變量X和Y的線性相關(guān)的程度，并常用r來表示，r的值介于-1，1之間，它的意義由圖10可以知道。r的絕對值越接近于1表示X和Y之間的線性關(guān)系越密切；r0，兩者呈正比關(guān)系，叫正相關(guān)；r 0兩者呈負相關(guān)。r的值接近于0，兩者沒有線性相關(guān)關(guān)系。圖10中(c)表示X和Y沒有任何關(guān)系，(d)表示X和Y有非線性相關(guān)關(guān)系，r的計算公式為 (2.6)式中 (2.7)對例3 =764.861r=415.605/ =0.9894 后者很接近于1，故最大積雪深度與灌溉面積有很密切的線性相關(guān)關(guān)系，且是正相關(guān).但是，相關(guān)系數(shù)有一個缺點，就是它接近1的程度與樣本的組數(shù)n是有關(guān)的，當n較小時，相關(guān)系數(shù)的絕對值容易接

6、近于1，當n較大時，相關(guān)系數(shù)的絕對值容易偏小。特別當n=2時，因為兩點決定一條直線，所以相關(guān)系數(shù)的絕對值總為1，在許多統(tǒng)計書中29給出相關(guān)系數(shù)的起碼值，當相關(guān)系數(shù)的絕對值大于表中之值時才可以認為X和Y有線性關(guān)系。此例當顯著性水平=1%時，表中的起碼值為0.765，今計算r=0.9894 0.765，故最大積雪深度與灌溉面積有高度的線性關(guān)系。在有些統(tǒng)計軟件中，常給出，這時便于區(qū)別記為。 (b)方差分析和F檢驗 因變量的波動可用來表達，這個波動是由兩個因素造成的；一個是X的變化引起Y相應(yīng)的變化，另一個是隨機誤差。前者造成Y的波動可用回歸平方和來表達，后者用殘差平方和來度量。它們分別用 和來表示，從

7、數(shù)學(xué)上可以導(dǎo)出 - (2.8)當X和Y為線性回歸模型(2.1)時，它們有如下更方便的計算公式- (2.9) 利用統(tǒng)計量F (2.10)可以來檢驗回歸方程(2.1)是否可信.當方程可信時F ，這里為F表中的臨界值，1和n-2為自由度，為顯著水平.對例3可以算得=1.802×415.606=748.922=764.961-748.922=16.039 F=8×748.922/16.039=373.55當=1%時。用F值和F表上的臨界值相比，若F ，表明Y的變化主要是由X的變化造成的，回歸方程(2.1)可信；若F值小于，回歸方程不可信。可信的程度也可分成不同等級，在本書中，=5%

8、時可信用“*” 表示，=1%時可信用“*” 表示。上述計算結(jié)果常列成方差分析表，如表10所示。表10 方差分析表 方差來源平方和自由度均方F顯著性回歸748.9221748.922373.550*誤差16.03982.005總和764.9619 (c) 殘差分析稱為殘差，它能提供許多有用的信息，表11給出了例3的10個殘差，利用殘差可以提供如下信息： (i)之估計 (2.11)給出了回歸方程的精度，它稱為殘差標準差，若隨機誤差遵從正態(tài)分布N(0，)，則Y的預(yù)報落在之內(nèi)的概率大約為95%，對例3可以算得=1.416，且10個均落于2×1.416之內(nèi)。 (ii)數(shù)據(jù)和模型之診斷 由殘差之

9、大小，可以發(fā)現(xiàn)異常(或叫離群)數(shù)據(jù)，可以發(fā)現(xiàn)模型(2.1)是否合適，是否要用非線性回歸模型等，這些已形成一整套理論，稱為回歸診斷，有興趣的讀者可參見文獻31。2.2多元線性回歸模型 當影響因變量Y的自變量不止一個時，比如有m個，這時Y和X之間的線性回歸方程為 (2.12)其中為回歸系數(shù)，為隨機誤差，常假定 。 設(shè)為觀測值，回歸分析的首要任務(wù)是利用它們來估計和，它們的最小二乘估計記作求估計值 需要解下面的線性方程組 (2.13)其中 當 求得后，計算 (2.14)回歸方程(2.12)建立后，檢驗其是否可信可用方差分析，這時公式(2.8)依然有效，但方差分析表(參看表10)將成為表12之形式，其中

10、 (2.15) 表12 方差分析表方差來源平方和自由度均方F顯著性回歸m殘差n-m-1總和n-1 它將與F的臨界值來比較，其比較的結(jié)果和結(jié)論請參見上節(jié)的討論，反映回歸精度的的估計公式為 (2.16) 類似于一元回歸相產(chǎn)系數(shù)r，可以定義適用于多元回歸的全關(guān)系數(shù)R，R定義)為和的相關(guān)系數(shù)，或定義為 (2.17) 例4 試用線性回歸模型(2.10)來擬合表9的試驗數(shù)據(jù)。 解：這時n=7，7組觀察值為(0.330，1.0，13，1.5)，(0.336，1.4，19，3.0)， (0.482，3.4，28，3.5)，它們的均值和為由于，故它們不必全部列出，將它們代入到方程級(2.13)中可以解得從而&#

11、215;2.2+0.00343 ×19-0.077×2.0 =0.201的估計為.于是回歸方程為 (2.18)進一步對它作方差分析，其方差分析表列于表13.表13 方差分析表方差來源自由度平方和均方F回歸30.0487700.0162573.29誤差30.0148380.004946總和60.063608當 =0.05 時F表的臨界值，回歸方程(2.18)不可信.這時，是否Y和三個因素之間不可能建立回歸關(guān)系呢？不是的，我們還應(yīng)作進一步探討，在下節(jié)我們將繼續(xù)討論該例。2.3 二次型回歸模型與變量篩選 由于因變量常常有交互作用，回歸模型(2.12)不足以反映實際，于是二次型回歸

12、模型常常為人們所采用.若有m個因素則二次型回歸模型為 (2.19)其中為回歸系數(shù)，為隨機誤差.我們看到，這時除了常數(shù)項 以外，方程有m(m+3)/2 項，當m=1，2， 時項數(shù)為m12345678910項數(shù)25914202735445465若使回歸系數(shù)的估計有可能，必要條件為n>1+m(m+3)/2.當m 較大時，通常不能滿足這個必要條件.于是有必要從方程(2.19)中選擇貢獻顯著的項，刪除不重要的項。有時，實際問題需要考慮高階的交互作用，如 等，這時篩選變量的任務(wù)就更為重要。在回歸分析中，有許多有效的篩選變量的技術(shù)，如a) 前進法，b) 后退法，c) 逐步回歸法，d) 最優(yōu)子集法(參看

13、25)。本章僅僅采用逐步回歸技術(shù)來篩選變量，這并不意味著逐步因歸在上述四項技術(shù)中最好的。 逐步回歸是回歸分析中的一種篩選變量的技術(shù)。開始它將貢獻最大的一個變量選入回歸方程，并且預(yù)先確定兩個閾值 和 ，用于決定變量能否入選或剔除。逐步回歸在每一步有三種可能的功能： a) 將一個新變量引進回歸模型，這時相應(yīng)的F統(tǒng)計量必須大于 b) 將一個變量從回歸模型中剔除，這時相應(yīng)的F統(tǒng)計量必須小于 c) 將回歸模型內(nèi)的一個變量和回歸模型外的一個變量交換位置。 執(zhí)行功能a)和b)時要注意如下原則： 設(shè)在當前步驟中有s個變量不在回歸模型中，有t個變量在回歸模型中。今欲從s個變量中挑選一個加入回歸模型之中，顯然應(yīng)挑

14、選使回歸效果最好的變量。這里回歸的效果可用方差分析表(見表10，表13)中F值來衡量，顯然我們要從s個變量中挑選一個變量使F值達到極大。類似地，若欲從t 個變量中刪除一個變量使其離開回歸模型，我們就是要選擇刪除后使回歸效果最好的變量，或選擇對當前回歸模型貢獻最小的變量。如果在某一步中，既能實現(xiàn)a)又能實現(xiàn)b)，兩者之和就是功能c)。大部分統(tǒng)計軟件包均有逐步回歸之功能，例如中國均勻設(shè)計學(xué)會推薦的軟件包。 現(xiàn)在我們對例4繼續(xù)進行討論。設(shè)先用后退法來選變量。所謂后退法，就是開始將所有的變量全部采用，然后逐步剔除對方程沒有顯著貢獻的變量，直到方程中所有的變量都有顯著貢獻為止。 仍考慮線性模型，開始三個

15、因素全部進入方程，得(2.18)。統(tǒng)計軟件包通常還會提供每個變量的t值，t值越大(按絕對值計)表示該因素越重要。對例2有這表明三個因素中以(反應(yīng)時間)對得率(Y)影響最大，配比次之，吡啶量最小。這些t 值都是隨機變量，它們遵從分布。若取=0.05 ，這時n=7，m=3， = 的臨界值(0.05)=3.18。t 值大于該值的因素表示對方程有顯著貢獻，否則表示不顯著。今 均小于(0.05)=3.18 ，說明回歸方程(2.18)的三個變量至少有一個不起顯著作用。于是我們將貢獻最小的刪去，重新建立Y和及的線性回歸方程，得 (2.20)，三個t 值分別為這時這三個t值遵從含四個自由度的t 分布，臨界值為

16、(0.05)=2.78，從而 應(yīng)從方程中剔除。然后對Y和建立回歸方程 (2.21)相應(yīng)的。因此，回歸方程(2.21)為“最終”的回歸模型。這里最終加上引號，表示并非真正的最終模型，而是在線性模型框架下的最終產(chǎn)物。 上述的分析只發(fā)現(xiàn)對Y有顯著作用，其它兩個因素均沒有顯著作用，該結(jié)論與實際經(jīng)驗不吻合，因此，猜想用線性模型不一定符合實際。于是進一步考試二次回歸模型(2.19)。這時方程中有9項(不算)。利用逐步回歸技術(shù)求得回歸方程如下： (2.22)其相應(yīng)的 。顯然，回歸方程(2.22)的效果優(yōu)于回歸方程(2.21)。方程(2.22)表明，因素和交互作用對Y有顯著的影響。值得注意的是，有些人對回歸分

17、析沒有足夠的理解，片面追求大的(或小的)，致使選進方程中的項過多，使誤差自由度為1或甚至為0，這時有關(guān)的結(jié)可靠性是很差的。因此，不應(yīng)片面追求大的，應(yīng)選擇n 稍大的均勻設(shè)計表，使得誤差有足夠的自由度5。2.4 應(yīng)用實例 均勻設(shè)計和正交設(shè)計以及其他試驗設(shè)計方法一樣，在工農(nóng)業(yè)生產(chǎn)和科學(xué)實驗中有廣闊的應(yīng)用前景，本文的文獻中列舉了部分應(yīng)用成果，其中有些成果成績顯著。更多的成果有待搜集。均勻設(shè)計不但在國內(nèi)得到廣泛應(yīng)用，在彼得格勒和香港也已開花結(jié)果。本節(jié)選擇香港浸會學(xué)院生物系的一項試驗，供讀者參考。例5 為了研究環(huán)境污染對人體的危害，今考核六種金屬的含量：鎘(Cd)，銅(Cu)，鋅(Zn)，鎳(Ni)，鉻(

18、Cr)，鉛(Pb)，每種金屬含量分別取了17個水平(百萬分之一，ppm)：0.01，0.05，0.1，0.2，0.4，0.8，1，2，4，5，8，10，12，14，16，18，20。今欲考慮這些金屬含量(包括它們的交互作用)對老鼠壽命的影響，該試驗考核老鼠身上某種細胞的死亡率。它們選用表，根據(jù)使用表的指示，它們選用了表中1，4，6，10，14，15列來安排六個因素，其試驗方案如表14所示。試驗的結(jié)果為死亡率。為了了解試驗誤差，提高結(jié)論的精度，他們在同一試驗條件下將試驗重復(fù)三次，三次結(jié)果()列于表15，三次死亡率的均值為，列于表15的最后一列。我們看到第17號試驗的死亡率為最高，因為這時六種金屬

19、都是最高含量，表明這些金屬對老鼠細胞確有致命作用。 現(xiàn)進一步用回歸分析來分析數(shù)據(jù)。由于每種金屬的含量由0.01(ppm)變到20(ppm)，最大得出小相差2000倍，于是直接用各因素的水平值作回歸不易獲得好的結(jié)果，通常要對水平值先作變換，用變換后的數(shù)據(jù)進行回歸。最常見的變換是取對數(shù)。于是回歸分析中的自變量成為logCd，logCu，logZn，logNi，logCr和logPb。根據(jù)以往經(jīng)驗，知道六種金屬間有交作用，故應(yīng)選用二次型回歸模型，并用逐步回歸來篩選變量。用同樣的 和 ，對 和分別進行逐步回歸，發(fā)現(xiàn)四組數(shù)據(jù)的結(jié)果非常吻合，表明試驗誤差不大，該試驗可以獲得可靠結(jié)論。為節(jié)省篇幅，我們僅列出

20、對的回歸方程=32.68+5.03LogCd+3.84LogCu+2.03LogNi+0.55(LogCu) -0.63(LogZn)+0.94(LogNi)+0.53(LogCd)(LogCu) -0.70(LogCd)(LogCr)+0.92(LogCu)(LogPb)方程中每一項的t 值分別為(常數(shù)項t 值未列)11.5，7.8，4.9，2.6，-3.4，4.1，2.4，-2.8，5.3，它們均遵從 分布，因(0.05)=2.365 小于上述所有t值之絕對值，故方程可信。表15 死亡率17.9517.6518.3317.922.0922.8522.6222.531.7432.7932.8

21、732.439.3740.6537.8739.331.9031.1833.7532.231.1430.6631.1831.039.8139.6140.8040.042.4841.8643.7942.724.9724.6525.0524.850.2951.2250.5450.660.7160.4359.6960.267.0171.9967.1268.732.7730.8633.7032.429.9428.6830.6629.767.8769.2567.0468.055.5655.2856.5255.779.5779.4378.4879.1 由方程我們可以給出如下結(jié)論：a)Cd，Cu 和Ni含量過

22、高，對老鼠細胞的死亡率有顯著作用，b)金屬Cd和Cu，Cd和Cr，Cu和Pb有交互作用，其中Cd和Cu，Cu和Pb對死亡率起正交互作用，而Cd和Cu對死亡率起負交互作用，c)Zn可能會中和其它金屬的破壞作用，降低老鼠細胞的死亡率。2.5 尋求最優(yōu)工藝條件 試驗設(shè)計的目的通常主要有二個，一是揭示變量(Y)與各因素之間的定性關(guān)系，二是尋求最優(yōu)工藝條件，回歸方程的建立可以達到一箭雙雕的目的。 現(xiàn)以例2來說明如何尋求最好的工藝條件，表9告訴我們，第7號試驗是7次試驗中最好工藝條件，即配比3.4，吡啶量28，反應(yīng)時間3.5 ，這個工藝條件和最優(yōu)工藝條件常常是很接近的。 在上述討論中，我們最終建立了回歸模

23、型(2.22)。該方程一般僅在試驗范圍內(nèi)成立，即配比，吡啶量10-28，反應(yīng)時間。尋求最優(yōu)模型等價于在這個范圍內(nèi)求方程(2.22)中的極大值。如果回歸方程比較復(fù)雜，可以用任何一個優(yōu)化算法(參見文獻33，34)來求最佳工藝條件，許多軟件包都含有優(yōu)化算法。數(shù)論方法也可以用來求的極大值，方開泰和王元提出了一個序貫算法SNTO，可以方便地求得的極大值，鑒于篇幅，這里就不詳細介紹了，有興趣的讀者可以參看文獻16。 對例2來講，可以用簡單的微積分求得極值，由于X在試驗范圍內(nèi)恒正，故由(2.22)知X 越大，越高，故X應(yīng)取試驗范圍內(nèi)極大值3.4。將X=3.4 代入(2.22)得令，解得=0，=2.7575，這時的極大值為51.85%。工藝條件=3.4，=2.7575 并未出現(xiàn)在原有試驗方案中，故應(yīng)在這個條件追加試驗，由于的最佳條件在試驗范圍邊界，故應(yīng)擴大試驗范圍。 對于許多實際工作者，不一定熟悉優(yōu)化方法，手邊沒有優(yōu)化的軟件。他們也不一定知道SNTO，也可能不會用微積分去求解極值。下面介紹一種“笨”辦法，其計算量較大，但程序好編。這種方法將每個因素的試驗