加性與廣義加性模型回歸分析

胡良平

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

1 概 述［1-2］

1.1 加性模型

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

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

1.2 廣義加性模型

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

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

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

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

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

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

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

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

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

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

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

1.4.2 后退擬合算法

1.4.2.1 未加權的后退擬合算法

未加權的后退擬合算法步驟如下：

第1步：初始化。

第2步：迭代。令 m=m+1；讓 j從1到 p循環，循環體內的計算公式為：

第3步：終止。

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

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

1.4.2.2 加權的后退擬合算法

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

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

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

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

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

表1 局部計分算法涉及到的重要內容

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

1.4.4 一般的局部計分算法

第1步：初始化。

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

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

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

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

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

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

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

2.1 問題與數據結構

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

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

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

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

2.2 基于常規方法構建多重線性回歸模型［3-4］

2.2.1 創建SAS數據集

創建一個名為“ExperimentA”的臨時SAS數據集所需的數據步程序：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程序說明】數據中每行上有3次試驗的結果，每次試驗結果都有3個數據，即溫度數值（temperature）、催化劑量（catalyst）與產量（yield）。

創建一個名為“ExperimentB”的臨時SAS數據集的SAS數據步程序：

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程序產生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 基于常規方法構建多重線性回歸模型

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

proc reg data=ExperimentA；

model y=x1 x2／r；

run；

記以上SAS程序創建的二重線性回歸模型為模型（1）。

proc reg data=ExperimentA；

model y=x1 x2／noint r；

run；

記以上SAS程序創建的二重線性回歸模型為模型（2）。

經比較，模型（1）優于模型（2）。具體方法詳見下文，此處從略。

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

proc reg data=ExperimentB；

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

run；

記以上SAS程序創建的多重線性回歸模型為模型（3）。

proc reg data=ExperimentB；

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

run；

記以上SAS程序創建的多重線性回歸模型為模型（4）。

經比較，模型（4）優于模型（3）。具體方法詳見下文，此處從略。

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

proc reg data=ExperimentB；

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

run；

記以上SAS程序創建的多重線性回歸模型為模型（5）。

proc reg data=ExperimentB；

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

記以上SAS程序創建的多重線性回歸模型為模型（6）。

經比較，模型（5）優于模型（6）。具體方法詳見下文，此處從略。

將模型（4）與模型（1）比較，得出模型（4）優于模型（1）。最后，需要將模型（5）與模型（4）作比較，具體方法如下：

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

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

利用下面的F檢驗對上述回歸模型（5）與模型（4）進行擬合優度比較：

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

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

2.3 基于加性模型構建多重回歸模型［2］

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

proc gam data=ExperimentA；

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

output out=a3 residual；

run；

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

【SAS輸出結果及其解釋】

因變量：y

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

輸入數據集的匯總

觀測數 112

缺失觀測數 0

分布 Gaussian

關聯函數 Identity

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

迭代匯總和擬合統計量

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

回歸模型分析

以上給出的是“加性模型”中“參數分析部分”的結果，即

平滑模型分析

以上給出的是“加性模型”中“非參數分析部分”的結果，即

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

平滑模型分析

以上是關于加性模型中兩個非參數項（即樣條函數）的假設檢驗結果，兩項各占用了3個自由度，經卡方檢驗，說明兩個非參數項都具有統計學意義。

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

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

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

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

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

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

Obs F F3_103

1 4.81687 2.69284

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

【結論】本例以加性模型的回歸分析結果為優。