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

胡良平

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

1　概　　述

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

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

y=f(x1,x2)

(1)

(2)

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

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

一般來(lái)說(shuō)，反應(yīng)曲面回歸分析應(yīng)用在下列場(chǎng)合：

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

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

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

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

1.3.1.1等高線

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

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

1.3.1.2穩(wěn)定點(diǎn)

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

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

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

1.3.2計(jì)算原理

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

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

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

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

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

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

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

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

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

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

研究?jī)?nèi)容：以結(jié)果變量MBT為因變量，以反應(yīng)溫度(temp)和作用時(shí)間(time)為兩個(gè)試驗(yàn)因素(或自變量)，構(gòu)建“二次反應(yīng)曲面回歸模型”；通過(guò)此曲面模型，洞察其表現(xiàn)，即繪制出“等高線”圖、找出其“穩(wěn)定點(diǎn)”；若“穩(wěn)定點(diǎn)”為“極值點(diǎ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程序說(shuō)明】調(diào)用RSREG過(guò)程進(jìn)行“反應(yīng)曲面回歸分析”；“plots=(ridge surface)；”要求系統(tǒng)對(duì)反應(yīng)曲面進(jìn)行“嶺分析”，即呈現(xiàn)出“等高線圖”和“穩(wěn)定點(diǎn)”；“l(fā)ackfit”要求系統(tǒng)進(jìn)行“失擬檢驗(yàn)”，即檢驗(yà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)計(jì)算結(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)于整個(gè)二次反應(yīng)曲面回歸模型和其中各部分的假設(shè)檢驗(yàn)結(jié)果：總模型具有統(tǒng)計(jì)學(xué)意義(F=4.81，P=0.041)，線性部分(即Temp和Time的一次項(xiàng))具有統(tǒng)計(jì)學(xué)意義，兩個(gè)平方項(xiàng)和一個(gè)交叉乘積項(xiàng)均無(wú)統(tǒng)計(jì)學(xué)意義。

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

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

參數(shù)自由度估計(jì)值標(biāo)準(zhǔn)誤差t值Pr>|t|代碼數(shù)據(jù)的參數(shù)估計(jì)值截距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)曲面回歸模型中各項(xiàng)的假設(shè)檢驗(yàn)結(jié)果，各項(xiàng)對(duì)應(yīng)的P均>0.05。這個(gè)結(jié)果并不理想。

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

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

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

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

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

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

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

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

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

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

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