胡良平

（1.軍事科學院研究生院，北京 100850；2.世界中醫(yī)藥學會聯(lián)合會臨床科研統(tǒng)計學專業(yè)委員會，北京 100029 *通信作者：胡良平，E-mail：lphu812＠sina.com）

1 概 述［1-2］

1.1 加性模型

將多重線性回歸模型進行推廣，使其表達式成為下面的式（1）形式：

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

1.2 廣義加性模型

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

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

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

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

1.3 加性與廣義加性模型回歸分析應用的場合

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

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

1.4.1 加性模型回歸分析的計算原理

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

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

由式（5）可知：對于所有其他的項“j”（j≠k），在給定“｛s＾j（），j≠k｝”的估計值時，其觀測值提供了一種用于估計每個光滑函數(shù)“sk（）”的方法。依此做法求得結(jié)果的迭代過程被稱為“后退擬合算法”，此法最早由Friedman和Stuetzle提出。

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)的計算公式為：

第3步：終止。

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

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

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

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

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

此處，ε=10-8是缺省的界值。用戶可以通過修改模型語句中的選項“EPSILON=”來改變此界值。

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

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

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

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

1.4.4 一般的局部計分算法

第1步：初始化。

si=g［E（y）］，s01=s02=… =s0p=0，m=0

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

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

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

此處，ε=10-8是缺省的界值。用戶可以通過修改模型語句中的選項“EPSILON=”來改變此界值。

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

2 基于加性模型回歸分析解決實際問題［2］

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

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

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

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

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

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

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

創(chuàng)建一個名為“ExperimentA”的臨時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程序說明】數(shù)據(jù)中每行上有3次試驗的結(jié)果，每次試驗結(jié)果都有3個數(shù)據(jù)，即溫度數(shù)值（temperature）、催化劑量（catalyst）與產(chǎn)量（yield）。

創(chuàng)建一個名為“ExperimentB”的臨時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；

【SAS程序說明】以上SAS程序產(chǎn)生7個“派生變量”，它們分別為x1與x2兩個原始自變量的平方項、立方項、交叉乘積項，具體地說，x3=x21、x4=x22、x5=x1×x2、x6=x31、x7=x32、x8=x21×x2、x9=x22×x1。其中，由前面的SAS程序可知：x1=temperature、x2=Catalyst、y=Yield。

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

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

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）。具體方法詳見下文，此處從略。

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

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）。具體方法詳見下文，此處從略。

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

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）。具體方法詳見下文，此處從略。

將模型（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檢驗對上述回歸模型（5）與模型（4）進行擬合優(yōu)度比較：

對應的 F臨界值 F（（2，106）（0.01））＜4.82，因 F=34.222＞4.82，說明P＜0.01，故需要選擇參數(shù)多的回歸模型（5）。模型（5）的輸出結(jié)果如下：

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

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

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

proc gam data=ExperimentA；

model y=spline（x1）spline（x2）；

output out=a3 residual；

run；

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

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

因變量：y

平滑模型成分：spline（x1）spline（x2）

輸入數(shù)據(jù)集的匯總

觀測數(shù) 112

缺失觀測數(shù) 0

分布 Gaussian

關(guān)聯(lián)函數(shù) Identity

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

迭代匯總和擬合統(tǒng)計量

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

回歸模型分析

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

平滑模型分析

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

將模型（7）與模型（8）合并成一個模型，見模型（9）：

平滑模型分析

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

圖1左邊的曲線描述的是模型（8）中的第1項；圖1右邊的曲線描述的是模型（8）中的第2項。其中，第1項比第2項更復雜。

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

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

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

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

已知：模型（9）的 SSε=68.464856、dfε=103；模型（5）的 SSε=78.07028、dfε=106。利用下面的SAS程序可以求出檢驗統(tǒng)計量F的數(shù)值以及對應的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é)果】

Obs F F3_103

1 4.81687 2.69284

因 F=4.817＞F（3，103）（0.95）=2.693，所以，P＜0.05，說明不能用含參數(shù)個數(shù)少的模型（5）取代含參數(shù)個數(shù)多的模型（9）。

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