胡良平

(1.軍事科學(xué)院研究生院，北京　100850；2.世界中醫(yī)藥學(xué)會聯(lián)合會臨床科研統(tǒng)計學(xué)專業(yè)委員會，北京　100029

1　概　　述

1.1　何為反應(yīng)曲面回歸分析

就是在構(gòu)建二重線性回歸模型時，引入兩個自變量的平方項和交叉乘積項，所構(gòu)建的二重線性回歸模型實際上是三維空間(x1，x2，y)中的一個二次曲面。由于此曲面的縱軸代表的是因變量(常被稱為“響應(yīng)變量”)，故稱此曲面為“響應(yīng)曲面”或“反應(yīng)曲面”。也就是說，y是關(guān)于(x1，x2)的二次函數(shù)。若自變量的個數(shù)多于2個，所形成的二次函數(shù)就應(yīng)該被稱為“超反應(yīng)曲面回歸模型”了。二次反應(yīng)曲面回歸模型如下面的式(1)和(2)所示：

y=f(x1,x2)

(1)

(2)

在式(1)和(2)中，y為定量的結(jié)果變量，即“因變量”或“響應(yīng)變量”；x1、x2為兩個原因變量，即自變量，在實際的多因素試驗研究中，它們通常是對試驗結(jié)果y有重要影響的試驗條件。

1.2　反應(yīng)曲面回歸分析應(yīng)用的場合

一般來說，反應(yīng)曲面回歸分析應(yīng)用在下列場合：

在一個化學(xué)或物理或生物學(xué)試驗研究中，涉及到兩個或兩個以上定量影響因素(即自變量)，當(dāng)各定量因素分別取一個特定水平時就構(gòu)成了一個特定的“試驗條件”。設(shè)有一個定量評價指標(biāo)y，則在任何一個特定條件下做試驗，y就會有一個取值。假定研究者在n個特定試驗條件下分別做了m次獨立重復(fù)試驗(若m=1就代表各試驗條件下只做了一次試驗)，此時，研究者最希望得到的結(jié)果是：①用一個二次方程式來定量反映因變量y隨自變量變化而變化的依賴關(guān)系；②各自變量分別取什么值(或水平)時，因變量y可以取得“最大值”或“最小值”。以前述的兩點為分析目的的“回歸分析”被稱為“反應(yīng)曲面回歸分析”[1]。

簡而言之，反應(yīng)曲面回歸分析常應(yīng)用于確定一個定量多因素的試驗研究中的“最佳生產(chǎn)條件”或“最佳工藝配方”，即在各定量因素分別取什么水平時做試驗，其試驗結(jié)果的取值最好(高優(yōu)指標(biāo)時，希望結(jié)果取得最大值，如產(chǎn)量；低優(yōu)指標(biāo)時，希望結(jié)果取得最小值，如能源消耗量)。

1.3　反應(yīng)曲面回歸分析的原理

1.3.1相關(guān)基本概念

1.3.1.1等高線

顧名思義，等高線就是在距離某個水平面相等高度上繪出的一條線。可以設(shè)想：一個技術(shù)非常高超的飛行員在距離某島嶼一萬米的高空中繞著此島嶼飛行一周，飛機尾部噴出的白色霧所形成的“圖案”就是一條等高線(高度為一萬米)。

用數(shù)學(xué)語言可描述如下：在一個三維直角坐標(biāo)系中，把一個二次反應(yīng)曲面想像成一個“圓頂草帽”，若用一把鋒利的鋸子在距離坐標(biāo)平面k個單位的高度上去平行地切割“圓頂草帽”(即具有高低起伏的曲面)，其切口就是一個環(huán)形的曲線，它就是該曲面在y=k時的等高線；當(dāng)k在y的取值范圍內(nèi)取一系列數(shù)值時，就形成了一系列的等高線。研究者根據(jù)這些等高線的形狀，就可比較清楚地看出：反應(yīng)曲面在二維平面上呈現(xiàn)出來的且在各個方向上的變化情況。

1.3.1.2穩(wěn)定點

反應(yīng)曲面上的穩(wěn)定點包括極大值點、極小值點和鞍點。

何為極大值點？若(x10，x20)為極大值點，在包含(x10，x20)的任何一個小區(qū)域內(nèi)的任何一點(x1i，x2i)，都滿足f(x10，x20)>f(x1i，x2i)，這里，y=f(x10，x20)就是(x1，x2)=(x10，x20)時y能取到的極大值；若前面所提及的“小區(qū)域”是兩個自變量變化的整個區(qū)域，則此極大值就是最大值。

同理，若(x10，x20)為極小值點，在包含(x10，x20)的任何一個小區(qū)域內(nèi)的任何一點(x1i，x2i)，都滿足f(x10，x20)f(x1i，x2i)或f(x10，x20)

1.3.2計算原理

反應(yīng)曲面回歸模型的構(gòu)建，若采用SAS中REG過程來實現(xiàn)，需要引入全部自變量的二次項和交叉乘積項，參見文獻(xiàn)[2]；若采用SAS中RSREG過程來實現(xiàn)，則非常簡單(參見表1后面的SAS程序)。仍可以采用最小平方法原理推導(dǎo)出正規(guī)方程組，通過求解此方程組可以獲得模型(1)中參數(shù)(即截距項和回歸系數(shù))的估計值。

進(jìn)行反應(yīng)曲面回歸分析的關(guān)鍵是如何對所構(gòu)建的二次反應(yīng)曲面回歸模型進(jìn)行分析，它涉及到“等高線”“穩(wěn)定點”等的計算，因篇幅所限，需要時可查閱有關(guān)文獻(xiàn)[3]。

2　基于SAS進(jìn)行反應(yīng)曲面回歸分析

2.1　問題與數(shù)據(jù)結(jié)構(gòu)

設(shè)有一個化學(xué)試驗，涉及到反應(yīng)溫度(temp)和作用時間(time)，受試對象為用于化學(xué)試驗的“樣品”。當(dāng)溫度和時間分別取某特定值時，就構(gòu)成了一個特定的試驗條件，試驗之后，定量的試驗結(jié)果“巰基苯并噻唑”(MBT)就會有一個具體的取值。某研究者考慮了“4.0、6.3、12.0、17.7和20.0”(h) 5個不同的反應(yīng)時間，又考慮了“220、229、250、271和280”(℃)5個不同的溫度。兩個試驗因素全部水平組合共有25種，某研究者只選取了其中一部分試驗條件進(jìn)行試驗，其因素水平組合及其試驗結(jié)果見表1。

表1　在不同反應(yīng)時間和溫度條件下進(jìn)行某化學(xué)試驗 得到某種物質(zhì)的產(chǎn)率MBT的結(jié)果

注：表1對應(yīng)的試驗及資料摘自SAS 9.3軟件RSREG過程的第1個“樣例”

【說明】在表1的每一行中，反應(yīng)時間與溫度分別取不同數(shù)值時，構(gòu)成一個特定的試驗條件(被稱為試驗點或設(shè)計點)，完全不同的試驗條件(試驗點或設(shè)計點)只有9個，因為(4.0，250)出現(xiàn)了兩次、(12.0，250)出現(xiàn)了三次。

2.2　研究目的與內(nèi)容

研究目的：試通過分析表1中的資料，求出在反應(yīng)時間和溫度分別取什么數(shù)值條件下，所得到的產(chǎn)率MBT最高。這樣的研究目的常被稱為“最優(yōu)生產(chǎn)條件的確定”問題。

研究內(nèi)容：以結(jié)果變量MBT為因變量，以反應(yīng)溫度(temp)和作用時間(time)為兩個試驗因素(或自變量)，構(gòu)建“二次反應(yīng)曲面回歸模型”；通過此曲面模型，洞察其表現(xiàn)，即繪制出“等高線”圖、找出其“穩(wěn)定點”；若“穩(wěn)定點”為“極值點”，再進(jìn)一步求出具體的“極大值”或“極小值”。

2.3　所需要的SAS程序

所需要的SAS程序如下：

data a;

input Time Temp MBT;

label Time=“Reaction Time (Hours)”

Temp=“Temperature (Degrees Centigrade)”

MBT=“Percent Yield Mercaptobenzothiazole”;

datalines;

(此處輸入表1中12行3列數(shù)據(jù))

;

run;

ods graphics on;

proc rsreg data=a plots=(ridge surface);

model MBT=Time Temp / lackfit;

ridge max;

run;

ods graphics off;

【SAS程序說明】調(diào)用RSREG過程進(jìn)行“反應(yīng)曲面回歸分析”；“plots=(ridge surface)；”要求系統(tǒng)對反應(yīng)曲面進(jìn)行“嶺分析”，即呈現(xiàn)出“等高線圖”和“穩(wěn)定點”；“l(fā)ackfit”要求系統(tǒng)進(jìn)行“失擬檢驗”，即檢驗可否用“多重線性回歸模型”取代“二次反應(yīng)曲面回歸模型”；“ridge max”要求系統(tǒng)求出響應(yīng)曲面上因變量的“最大值”。

2.4　SAS輸出結(jié)果及其解釋

變量“MBT”的響應(yīng)曲面:PercentYield Mercaptobenzothiazole響應(yīng)均值79.916667均方根誤差4.615964R20.8003偏差系數(shù)5.7760

以上是關(guān)于因變量MBT的有關(guān)計算結(jié)果。

回歸自由度I型平方和R2F值Pr>F線性2313.5858030.48997.360.0243二次2146.7681440.22933.440.1009叉積151.8400000.08102.430.1698總模型5512.1939470.80034.810.0410

以上是關(guān)于整個二次反應(yīng)曲面回歸模型和其中各部分的假設(shè)檢驗結(jié)果：總模型具有統(tǒng)計學(xué)意義(F=4.81，P=0.041)，線性部分(即Temp和Time的一次項)具有統(tǒng)計學(xué)意義，兩個平方項和一個交叉乘積項均無統(tǒng)計學(xué)意義。

殘差自由度平方和均方F值Pr>F缺少擬合3124.69605341.56535139.630.0065純誤差33.1466671.048889總誤差6127.84272021.307120

以上是關(guān)于“失擬檢驗”的結(jié)果：線性不能描述的部分(即擬合失敗部分)具有統(tǒng)計學(xué)意義(F=39.63，P=0.0065)，說明有必要構(gòu)建二次反應(yīng)曲面回歸模型。

參數(shù)自由度估計值標(biāo)準(zhǔn)誤差t值Pr>|t|代碼數(shù)據(jù)的參數(shù)估計值截距1-545.867976277.145373-1.970.096482.173110Time16.8728635.0049281.370.2188-1.014287Temp14.9897432.1658392.300.0608-8.676768Time*Time10.0216310.0567840.380.71641.384394Temp*Time1-0.0300750.019281-1.560.1698-7.218045Temp*Temp1-0.0098360.004304-2.290.0623-8.852519

以上呈現(xiàn)出二次反應(yīng)曲面回歸模型中各項的假設(shè)檢驗結(jié)果，各項對應(yīng)的P均>0.05。這個結(jié)果并不理想。

因子自由度平方和均方F值Pr>F標(biāo)簽Time361.29095720.4303190.960.4704Reaction Time (Hours)Temp3461.250925153.7503087.220.0205Temperature (Degrees Centigrade

以上是對涉及“Time”的三項合并的總評價(無統(tǒng)計學(xué)意義，F(xiàn)=0.96，P=0.4704)、對涉及“Temp”的三項合并的總評價(F=7.22，P=0.0205)，此結(jié)果表明：僅溫度(Temp)對試驗結(jié)果(MBT)的影響具有統(tǒng)計學(xué)意義。

自變量不同水平組合下因變量MBT的最大值及其標(biāo)準(zhǔn)誤見表2。

表2　自變量不同水平組合下因變量MBT的最大值及其標(biāo)準(zhǔn)誤

由表2可知：當(dāng)Time = 18.451、Temp = 232.256時，因變量MBT = 87.733為各種水平組合條件下的“最大值”。

圖1　具有設(shè)計點的“MBT”的響應(yīng)等高線

注：圖1左邊縱坐標(biāo)軸為“Reaction Time(h)”，即“反應(yīng)時間(h)”；右邊縱坐標(biāo)軸為“標(biāo)準(zhǔn)誤差”；橫坐標(biāo)軸為“Temperature(Degrees Centigrade)”，即“溫度(攝氏度)”

在圖1中出現(xiàn)了9個“圓圈”，它們代表表1中9個不同的“設(shè)計點(或試驗點)”。圖中的“弧線”反映了某些“設(shè)計點”及其附近“未做試驗的點”對應(yīng)的結(jié)果變量(MBT)的等高線?，F(xiàn)以圖1左上角標(biāo)注“90”的那條“弧線”為例，說明“等高線”的含義：“90”代表“結(jié)果變量(MBT)的取值為90”之意，“弧線”代表橫坐標(biāo)軸上的“溫度”的變化范圍大約在(220，235)之間，而縱坐標(biāo)軸上的“反應(yīng)時間”的變化范圍大約在20℃以上。讀者可以嘗試去解釋標(biāo)注“80”的兩條等高線的含義，此處不再贅述。

【說明】在上面的輸出結(jié)果中有這樣一個結(jié)果：“兩個平方項和一個交叉乘積項均無統(tǒng)計學(xué)意義”，這個結(jié)果提示：本例的試驗研究中所涉及的兩個定量因素對定量結(jié)果的影響不存在“二次項”效應(yīng)，換句話說，基于本例資料，不需要采用反應(yīng)曲面回歸模型，而只需要采用“一元二重線性回歸模型(用幾何學(xué)來理解，就是一個二維平面，而不是一個三維曲面)”。若果真這樣做，就不存在“反應(yīng)曲面”了，也就不存在因變量的“最大值”或“最小值”了。由此得出如下推論：

(1)反應(yīng)曲面回歸分析最適合用于需要確定多定量因素的最佳生產(chǎn)條件的試驗研究場合。

(2)應(yīng)用此法的前提條件是：參與試驗研究的定量因素已經(jīng)過專業(yè)和統(tǒng)計學(xué)方法嚴(yán)格篩選并被保留下來，并且它們之間的全部或大多數(shù)二次項和交叉乘積項都具有統(tǒng)計學(xué)意義。