胡良平

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

1 概 述［1-2］

1.1 加性模型

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

在式（1）中，sj（Xj），j=1，2，…，p，是 P個(gè)“光滑函數(shù)”；誤差“ε”滿(mǎn)足如下條件：它的期望為0［E（ε）=0］、方差為 σ2［Var（ε）=σ2］。為了使式（1）成為可以估計(jì)的，要求光滑函數(shù) si（Xi）必須滿(mǎn)足如下的標(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）中，θ被稱(chēng)為“自然參數(shù)”，φ被稱(chēng)為“尺度參數(shù)”；以因變量y的平均值μ為自變量構(gòu)造出的函數(shù) g（μ）被稱(chēng)為“連接函數(shù)”，它與協(xié)變量“X1、X2、…、Xp”之間建立了關(guān)系。下面的數(shù)量定義了“加性分量（或成分）”，見(jiàn)式（3）：

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

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

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

廣義線(xiàn)性模型強(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）被稱(chēng)為“第k項(xiàng)光滑參數(shù)與因變量y”之間的“殘差”，即 Rk≈sk（Xk）。嚴(yán)格地說(shuō)，應(yīng)該有下式成立，見(jiàn)式（6）：

由式（5）可知：對(duì)于所有其他的項(xiàng)“j”（j≠k），在給定“｛s＾j（），j≠k｝”的估計(jì)值時(shí)，其觀(guān)測(cè)值提供了一種用于估計(jì)每個(gè)光滑函數(shù)“sk（）”的方法。依此做法求得結(jié)果的迭代過(guò)程被稱(chēng)為“后退擬合算法”，此法最早由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)的計(jì)算公式為：

第3步：終止。

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

值得注意的是：對(duì)于除正態(tài)分布之外的其他分布而言，具有權(quán)重的數(shù)值不穩(wěn)定可能會(huì)引起收斂問(wèn)題。即使當(dāng)算法收斂時(shí)，各個(gè)個(gè)體函數(shù)并不需要彼此完全不同，即便對(duì)于同一個(gè)擬合曲面來(lái)說(shuō)，由于協(xié)變量之間的依賴(lài)性會(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是缺省的界值。用戶(hù)可以通過(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可知：一旦分布被指定，相應(yīng)的“那些量”也就被定義了。于是，可按下面的步驟實(shí)施“局部計(jì)分算法”。

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

第1步：初始化。

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

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

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

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

此處，ε=10-8是缺省的界值。用戶(hù)可以通過(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é)果

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

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

2.2 基于常規(guī)方法構(gòu)建多重線(xiàn)性回歸模型［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；

【SAS程序說(shuō)明】以上SAS程序產(chǎn)生7個(gè)“派生變量”，它們分別為x1與x2兩個(gè)原始自變量的平方項(xiàng)、立方項(xiàng)、交叉乘積項(xiàng)，具體地說(shuō)，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)建多重線(xiàn)性回歸模型

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

proc reg data=ExperimentA；

model y=x1 x2／r；

run；

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

proc reg data=ExperimentA；

model y=x1 x2／noint r；

run；

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

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

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

proc reg data=ExperimentB；

model y=x1-x5／selection=backward sle=0.05 r；

run；

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

proc reg data=ExperimentB；

model y=x1-x5／noint selection=backward sle=0.05 r；

run；

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

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

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

proc reg data=ExperimentB；

model y=x1-x9／selection=backward sle=0.05 r；

run；

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

proc reg data=ExperimentB；

model y=x1-x9／noint selection=backward sle=0.05 r；run；

記以上SAS程序創(chuàng)建的多重線(xiàn)性回歸模型為模型（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é)果如下：

根據(jù)最后的“參數(shù)估計(jì)值”，請(qǐng)讀者寫(xiě)出相應(yīng)的“五重線(xiàn)性回歸模型”的表達(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ù)集的匯總

觀(guān)測(cè)數(shù) 112

缺失觀(guān)測(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ì)量

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

回歸模型分析

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

平滑模型分析

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

將模型（7）與模型（8）合并成一個(gè)模型，見(jiàn)模型（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左邊的曲線(xiàn)描述的是模型（8）中的第1項(xiàng)；圖1右邊的曲線(xiàn)描述的是模型（8）中的第2項(xiàng)。其中，第1項(xiàng)比第2項(xiàng)更復(fù)雜。

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

2.4 兩類(lèi)回歸模型擬合效果比較

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

將常規(guī)多重線(xiàn)性回歸模型中擬合效果最好的模型（5）與加性模型（9）進(jìn)行比較，用類(lèi)似于模型（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é)果】

Obs F F3_103

1 4.81687 2.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)。