胡良平

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

1 概 述[1-2]

1.1 加性模型

將多重線性回歸模型進(jìn)行推廣，使其表達(dá)式成為下面的式(1)形式：

Y=s0+s1(X1)+s2(X2)+…+sp(Xp)+ε

(1)

在式(1)中，sj(Xj),j=1,2,…,p,是P個(gè)“光滑函數(shù)”；誤差“ε”滿足如下條件：它的期望為0[E(ε)=0]、方差為σ2[Var(ε)=σ2]。為了使式(1)成為可以估計(jì)的，要求光滑函數(shù)si(Xi)必須滿足如下的標(biāo)準(zhǔn)化條件：即期望為0，E[sj(Xj)]=0。式(1)中的P個(gè)光滑函數(shù)不以參數(shù)形式呈現(xiàn)，而以非參數(shù)形式呈現(xiàn)。

1.2 廣義加性模型

在加性模型式(1)中，假定因變量y服從正態(tài)分布。然而，在很多場(chǎng)合下，因變量不服從正態(tài)分布，而可能服從其他某種分布?，F(xiàn)假定式(1)中的因變量y具有下面的指數(shù)族分布密度，見(jiàn)式(2)：

(2)

在式(2)中，θ被稱為“自然參數(shù)”，φ被稱為“尺度參數(shù)”；以因變量y的平均值μ為自變量構(gòu)造出的函數(shù)g(μ)被稱為“連接函數(shù)”，它與協(xié)變量“X1、X2、…、Xp”之間建立了關(guān)系。下面的數(shù)量定義了“加性分量(或成分)”，見(jiàn)式(3)：

(3)

式(3)中，S1()、S2()、…、Sp()都是“光滑函數(shù)”。μ與η之間的關(guān)系由下式來(lái)定義：

g(μ)=η

(4)

最常用的連接函數(shù)為“典型連接函數(shù)”，即η=θ。

1.3 加性與廣義加性模型回歸分析應(yīng)用的場(chǎng)合

廣義線性模型強(qiáng)調(diào)對(duì)模型中參數(shù)的估計(jì)和推斷，而廣義加性模型則聚焦于如何用非參數(shù)法探測(cè)數(shù)據(jù)。換言之，廣義加性模型更適合于探查數(shù)據(jù)并可視化因變量與自變量之間的關(guān)系。

1.4 加性與廣義加性模型回歸分析的計(jì)算原理

1.4.1 加性模型回歸分析的計(jì)算原理

基于加性模型式(1)，可以構(gòu)造如下的殘差，見(jiàn)式(5)：

(5)

式(5)被稱為“第k項(xiàng)光滑參數(shù)與因變量y”之間的“殘差”，即Rk≈sk(Xk)。嚴(yán)格地說(shuō)，應(yīng)該有下式成立，見(jiàn)式(6)：

E(Rk|Xk)=sk(Xk)

(6)

1.4.2 后退擬合算法

1.4.2.1 未加權(quán)的后退擬合算法

未加權(quán)的后退擬合算法步驟如下：

第1步：初始化。

第2步：迭代。令m=m+1；讓j從1到p循環(huán)，循環(huán)體內(nèi)的計(jì)算公式為：

第3步：終止。

直到上式的計(jì)算結(jié)果不再下降，或滿足收斂的臨界值，就停止迭代。

值得注意的是：對(duì)于除正態(tài)分布之外的其他分布而言，具有權(quán)重的數(shù)值不穩(wěn)定可能會(huì)引起收斂問(wèn)題。即使當(dāng)算法收斂時(shí)，各個(gè)個(gè)體函數(shù)并不需要彼此完全不同，即便對(duì)于同一個(gè)擬合曲面來(lái)說(shuō)，由于協(xié)變量之間的依賴性會(huì)導(dǎo)致多于一個(gè)表達(dá)式出現(xiàn)。

1.4.2.2 加權(quán)的后退擬合算法

除了要對(duì)光滑器進(jìn)行加權(quán)之外，加權(quán)的后退擬合算法具有與未加權(quán)的后退擬合算法相同的形式。在SAS的GAM過(guò)程步中，具體是在采用“局部計(jì)分過(guò)程”中且資料為非正態(tài)分布時(shí)，使用了“權(quán)重”。

GAM過(guò)程使用下面的“條件”作為后退擬合算法的“收斂臨界值”：

此處，ε=10-8是缺省的界值。用戶可以通過(guò)修改模型語(yǔ)句中的選項(xiàng)“EPSILON=”來(lái)改變此界值。

1.4.3 廣義加性模型回歸分析的計(jì)算原理

廣義加性模型回歸分析的計(jì)算原理比前面介紹的加性模型回歸分析的計(jì)算原理稍復(fù)雜，其中，最關(guān)鍵的內(nèi)容為“局部計(jì)分算法”。該算法的重要內(nèi)容取決于與每個(gè)特定分布對(duì)應(yīng)的“連接函數(shù)”。它們之間的關(guān)系見(jiàn)表1。

表1 局部計(jì)分算法涉及到的重要內(nèi)容

注：表1摘錄自SAS軟件GAM過(guò)程幫助信息

由表1可知：一旦分布被指定，相應(yīng)的“那些量”也就被定義了。于是，可按下面的步驟實(shí)施“局部計(jì)分算法”。

1.4.4 一般的局部計(jì)分算法

第1步：初始化。

第2步：迭代。令m=m+1；從前一次迭代中獲得各變量的相應(yīng)數(shù)值，這些變量分別是：預(yù)測(cè)量η、均值μ、權(quán)重w和校正后的因變量z：

第3步：終止。直到達(dá)到收斂臨界值或離差不再減少時(shí)，停止迭代。這里所說(shuō)的“離差”，實(shí)際上是廣義線性模型中“RSS”的一個(gè)擴(kuò)展或推廣。

GAM過(guò)程使用下面的“條件”作為局部計(jì)分算法的“收斂臨界值”：

此處，ε=10-8是缺省的界值。用戶可以通過(guò)修改模型語(yǔ)句中的選項(xiàng)“EPSILON=”來(lái)改變此界值。

算法小結(jié)：廣義加性模型的估計(jì)過(guò)程由兩個(gè)循環(huán)構(gòu)成。在局部計(jì)分算法(外循環(huán))每一步內(nèi)部，使用加權(quán)后退擬合算法(內(nèi)循環(huán))，直到收斂或RSS不再減少。然后，基于來(lái)自這個(gè)加權(quán)后退擬合算法得到的估計(jì)量，計(jì)算出一組新的權(quán)重，開(kāi)始計(jì)分算法的下輪迭代。當(dāng)達(dá)到收斂臨界值或估計(jì)量的離差停止下降時(shí)，計(jì)分算法也就停止了。

2 基于加性模型回歸分析解決實(shí)際問(wèn)題[2]

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

【例1】下面是一個(gè)假設(shè)的來(lái)自化學(xué)試驗(yàn)的例子：每次試驗(yàn)，研究者將某種催化劑加入到某種化學(xué)溶液中，從而合成一種新化合物。其數(shù)據(jù)是測(cè)量溶液的溫度(temperature)、加入的催化劑量(catalyst)和化學(xué)反應(yīng)量(yield)的結(jié)果。試驗(yàn)數(shù)據(jù)的結(jié)構(gòu)很簡(jiǎn)單，兩個(gè)計(jì)量原因變量和一個(gè)計(jì)量結(jié)果變量及其取值，前6次試驗(yàn)數(shù)據(jù)見(jiàn)表2。

表2 不同“溶液溫度”和“催化劑量”條件下化學(xué)反應(yīng)量的測(cè)定結(jié)果

注：表中僅列出了少量數(shù)據(jù)，詳細(xì)試驗(yàn)數(shù)據(jù)見(jiàn)下面的SAS程序，此處從略

【對(duì)數(shù)據(jù)結(jié)構(gòu)的分析】對(duì)于每次試驗(yàn)而言，可以觀測(cè)到3個(gè)計(jì)量的數(shù)據(jù)，即溶液的溫度(temperature)、加入的催化劑量(catalyst)和化學(xué)反應(yīng)的量(yield)。

【統(tǒng)計(jì)分析方法的選擇】若希望考察化學(xué)反應(yīng)的量(yield)是如何隨溶液的溫度(temperature)和催化劑量(catalyst)變化而變化的依賴關(guān)系，可選擇多重線性回歸分析。因本例中的因變量為計(jì)量變量，故可以考慮選用“加性模型回歸分析”。若因變量是定性變量或計(jì)數(shù)變量，就可能需要選用“廣義加性模型回歸分析”。因篇幅所限，本文只介紹如何用SAS實(shí)現(xiàn)“加性模型回歸分析”。

2.2 基于常規(guī)方法構(gòu)建多重線性回歸模型[3-4]

2.2.1 創(chuàng)建SAS數(shù)據(jù)集

創(chuàng)建一個(gè)名為“ExperimentA”的臨時(shí)SAS數(shù)據(jù)集所需的數(shù)據(jù)步程序：

data ExperimentA;

format Temperature f4.0 Catalyst f6.3 Yield f8.3;

input Temperature Catalyst Yield @@;

x1=temperature;x2=Catalyst;y=Yield;

datalines;

80 0.005 6.039 80 0.010 4.719 80 0.015 6.301

80 0.020 4.558 80 0.025 5.917 80 0.030 4.365

80 0.035 6.540 80 0.040 5.063 80 0.045 4.668

80 0.050 7.641 80 0.055 6.736 80 0.060 7.255

80 0.065 5.515 80 0.070 5.260 80 0.075 4.813

80 0.080 4.465 90 0.005 4.540 90 0.010 3.553

90 0.015 5.611 90 0.020 4.586 90 0.025 6.503

90 0.030 4.671 90 0.035 4.919 90 0.040 6.536

90 0.045 4.799 90 0.050 6.002 90 0.055 6.988

90 0.060 6.206 90 0.065 5.193 90 0.070 5.783

90 0.075 6.482 90 0.080 5.222 100 0.005 5.042

100 0.010 5.551 100 0.015 4.804 100 0.020 5.313

100 0.025 4.957 100 0.030 6.177 100 0.035 5.433

100 0.040 6.139 100 0.045 6.217 100 0.050 6.498

100 0.055 7.037 100 0.060 5.589 100 0.065 5.593

100 0.070 7.438 100 0.075 4.794 100 0.080 3.692

110 0.005 6.005 110 0.010 5.493 110 0.015 5.107

110 0.020 5.511 110 0.025 5.692 110 0.030 5.969

110 0.035 6.244 110 0.040 7.364 110 0.045 6.412

110 0.050 6.928 110 0.055 6.814 110 0.060 8.071

110 0.065 6.038 110 0.070 6.295 110 0.075 4.308

110 0.080 7.020 120 0.005 5.409 120 0.010 7.009

120 0.015 6.160 120 0.020 7.408 120 0.025 7.123

120 0.030 7.009 120 0.035 7.708 120 0.040 5.278

120 0.045 8.111 120 0.050 8.547 120 0.055 8.279

120 0.060 8.736 120 0.065 6.988 120 0.070 6.283

120 0.075 7.367 120 0.080 6.579 130 0.005 7.629

130 0.010 7.171 130 0.015 5.997 130 0.020 6.587

130 0.025 7.335 130 0.030 7.209 130 0.035 8.259

130 0.040 6.530 130 0.045 8.400 130 0.050 7.218

130 0.055 9.167 130 0.060 9.082 130 0.065 7.680

130 0.070 7.139 130 0.075 7.275 130 0.080 7.544

140 0.005 4.860 140 0.010 5.932 140 0.015 3.685

140 0.020 5.581 140 0.025 4.935 140 0.030 5.197

140 0.035 5.559 140 0.040 4.836 140 0.045 5.795

140 0.050 5.524 140 0.055 7.736 140 0.060 5.628

140 0.065 6.644 140 0.070 3.785 140 0.075 4.853

140 0.080 6.006

;

run;

【SAS程序說(shuō)明】數(shù)據(jù)中每行上有3次試驗(yàn)的結(jié)果，每次試驗(yàn)結(jié)果都有3個(gè)數(shù)據(jù)，即溫度數(shù)值(temperature)、催化劑量(catalyst)與產(chǎn)量(yield)。

創(chuàng)建一個(gè)名為“ExperimentB”的臨時(shí)SAS數(shù)據(jù)集的SAS數(shù)據(jù)步程序：

data ExperimentB;

set ExperimentA;

x3=x1*x1;x4=x2*x2;x5=x1*x2;x6=x3*x1;

x7=x4*x2;x8=x3*x2;x9=x4*x1;

run;

2.2.2 基于常規(guī)方法構(gòu)建多重線性回歸模型

利用下面的兩個(gè)SAS過(guò)程步程序可以創(chuàng)建兩個(gè)二重線性回歸模型：

proc reg data=ExperimentA;

model y=x1 x2/r;

run;

記以上SAS程序創(chuàng)建的二重線性回歸模型為模型(1)。

proc reg data=ExperimentA;

model y=x1 x2/noint r;

run;

記以上SAS程序創(chuàng)建的二重線性回歸模型為模型(2)。

經(jīng)比較，模型(1)優(yōu)于模型(2)。具體方法詳見(jiàn)下文，此處從略。

引入自變量的“二次項(xiàng)”，利用下面的兩個(gè)SAS過(guò)程步程序可以創(chuàng)建兩個(gè)多重線性回歸模型：

proc reg data=ExperimentB;

model y=x1-x5/selection=backward sle=0.05 r;

run;

記以上SAS程序創(chuàng)建的多重線性回歸模型為模型(3)。

proc reg data=ExperimentB;

model y=x1-x5/noint selection=backward sle=0.05 r;

run;

記以上SAS程序創(chuàng)建的多重線性回歸模型為模型(4)。

經(jīng)比較，模型(4)優(yōu)于模型(3)。具體方法詳見(jiàn)下文，此處從略。

引入自變量的“三次項(xiàng)”，利用下面的兩個(gè)SAS過(guò)程步程序可以創(chuàng)建兩個(gè)多重線性回歸模型：

proc reg data=ExperimentB;

model y=x1-x9/selection=backward sle=0.05 r;

run;

記以上SAS程序創(chuàng)建的多重線性回歸模型為模型(5)。

proc reg data=ExperimentB;

model y=x1-x9/noint selection=backward sle=0.05 r;

run;

記以上SAS程序創(chuàng)建的多重線性回歸模型為模型(6)。

經(jīng)比較，模型(5)優(yōu)于模型(6)。具體方法詳見(jiàn)下文，此處從略。

將模型(4)與模型(1)比較，得出模型(4)優(yōu)于模型(1)。最后，需要將模型(5)與模型(4)作比較，具體方法如下：

模型(4)的有關(guān)信息為：SSε=128.48055(模型誤差的離均差平方和)、dfε=108(誤差的自由度)；

模型(5)的有關(guān)信息為：SSε=78.07028(模型誤差的離均差平方和)、dfε=106(誤差的自由度)。

利用下面的F檢驗(yàn)對(duì)上述回歸模型(5)與模型(4)進(jìn)行擬合優(yōu)度比較：

對(duì)應(yīng)的F臨界值F((2,106)(0.01))<4.82，因F=34.222>4.82，說(shuō)明P<0.01，故需要選擇參數(shù)多的回歸模型(5)。模型(5)的輸出結(jié)果如下：

方差分析源自由度平方和均方F值Pr > F模型588.6661117.7332224.08<.0001誤差10678.070280.73651校正合計(jì)111166.73639

均方根誤差0.85820R方0.5318因變量均值6.14077調(diào)整 R方0.5097變異系數(shù)13.97550

參數(shù)估計(jì)值變量自由度參數(shù)估計(jì)值標(biāo)準(zhǔn)誤差t 值Pr > |t|Intercept1122.1032018.516226.59<.0001x11-3.459610.52226-6.62<.0001x310.033220.004826.89<.0001x411437.83096246.932355.82<.0001x61-0.000103550.00001460-7.09<.0001x71-17556.000003130.36119-5.61<.0001

根據(jù)最后的“參數(shù)估計(jì)值”，請(qǐng)讀者寫出相應(yīng)的“五重線性回歸模型”的表達(dá)式，此處從略。

2.3 基于加性模型構(gòu)建多重回歸模型[2]

利用下面的SAS程序可基于加性模型構(gòu)建多重回歸模型：

proc gam data=ExperimentA;

model y=spline(x1) spline(x2);

output out=a3 residual;

run;

【SAS程序說(shuō)明】以上SAS程序調(diào)用GAM過(guò)程擬合加性模型。模型語(yǔ)句等號(hào)右邊的兩項(xiàng)分別用“三次樣條函數(shù)”擬合自變量x1與x2。

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

因變量: y平滑模型成分: spline(x1) spline(x2) 輸入數(shù)據(jù)集的匯總觀測(cè)數(shù)112缺失觀測(cè)數(shù)0分布Gaussian關(guān)聯(lián)函數(shù)Identity

以上是關(guān)于數(shù)據(jù)集一般情況的描述，并告知：假定因變量y服從正態(tài)分布(或高斯分布)、采用恒等的關(guān)聯(lián)函數(shù)，實(shí)際上，就是沒(méi)有對(duì)y作任何變量變換。

迭代匯總和擬合統(tǒng)計(jì)量最終 Backfitting 迭代次數(shù)2最終 Backfitting 準(zhǔn)則2.673726E-14最終估計(jì)的偏差68.464845603Backfitting 算法收斂

以上是關(guān)于“迭代匯總和擬合統(tǒng)計(jì)量”的信息，關(guān)鍵是倒數(shù)第二行：最終估計(jì)的偏差為68.464846，此值相當(dāng)于通?；貧w分析給出的“模型誤差的離差平方法和”。

回歸模型分析參數(shù)估計(jì)值參數(shù)參數(shù)估計(jì)值標(biāo)準(zhǔn)誤差t 值Pr > |t|Intercept3.776180.453488.33<.0001Linear(x1)0.017650.003854.58<.0001Linear(x2)9.950913.342392.980.0036

以上給出的是“加性模型”中“參數(shù)分析部分”的結(jié)果，即

模型(7)

以上給出的是“加性模型”中“非參數(shù)分析部分”的結(jié)果，即

模型(8)

將模型(7)與模型(8)合并成一個(gè)模型，見(jiàn)模型(9)：

模型(9)

以上是關(guān)于加性模型中兩個(gè)非參數(shù)項(xiàng)(即樣條函數(shù))的假設(shè)檢驗(yàn)結(jié)果，兩項(xiàng)各占用了3個(gè)自由度，經(jīng)卡方檢驗(yàn)，說(shuō)明兩個(gè)非參數(shù)項(xiàng)都具有統(tǒng)計(jì)學(xué)意義。

圖1左邊的曲線描述的是模型(8)中的第1項(xiàng)；圖1右邊的曲線描述的是模型(8)中的第2項(xiàng)。其中，第1項(xiàng)比第2項(xiàng)更復(fù)雜。

圖1 兩個(gè)非參數(shù)項(xiàng)分別與x1、x2之間的函數(shù)曲線

2.4 兩類回歸模型擬合效果比較

常規(guī)多重線性回歸模型與加性多重回歸模型對(duì)同一個(gè)資料究竟誰(shuí)的擬合效果更好？這個(gè)問(wèn)題尚無(wú)公認(rèn)的評(píng)判方法，但可以近似地采用下面的方法進(jìn)行比較：

將常規(guī)多重線性回歸模型中擬合效果最好的模型(5)與加性模型(9)進(jìn)行比較，用類似于模型(5)與模型(4)比較的F檢驗(yàn)：

已知：模型(9)的SSε=68.464856、dfε=103；模型(5)的SSε=78.07028、dfε=106。利用下面的SAS程序可以求出檢驗(yàn)統(tǒng)計(jì)量F的數(shù)值以及對(duì)應(yīng)的F臨界值：

data abc;

v1=(78.07028-68.464856)/(106-103);

v2=68.464856/103;

F=v1/v2;

F3_103=FINV(0.95,3,103);

proc print data=abc;

var F F3_103;

run;

【SAS輸出結(jié)果】

ObsFF3_10314.816872.69284

因F=4.817>F(3,103)(0.95)=2.693，所以，P<0.05，說(shuō)明不能用含參數(shù)個(gè)數(shù)少的模型(5)取代含參數(shù)個(gè)數(shù)多的模型(9)。

【結(jié)論】本例以加性模型的回歸分析結(jié)果為優(yōu)。