胡良平

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

1 概 述

1.1 如何直觀分析散布圖的表現(xiàn)

在單個自變量的回歸分析中，為了選擇到合適的回歸分析模型，最簡單且最有效的方法是先繪制(X，Y)的散布圖。通?？梢罁?jù)散布圖中所有散點的分布情況，結(jié)合初等函數(shù)的表達式及其圖象[1]，嘗試擬合幾種最可能的曲線類型(也包括直線)，并通過擬合優(yōu)度比較，最終確定一種最合適的曲線回歸模型。然而，如圖1中散點所呈現(xiàn)的“形態(tài)”，確實難以作出合理的判斷。所有散點幾乎在一個長方形區(qū)域內(nèi)“星羅密布”，說它們近似呈“均勻分布”似乎比較合理。

圖1 大氣壓值(pressure)隨年份值(year)變化而變化的散布圖

1.2 變量變換的必要性

對于圖1中散點的分布情況，若直接擬合直線回歸模型似乎是很合理的，但在統(tǒng)計學(xué)上又是毫無價值的(因為R2=0.0071)。所以，應(yīng)采用合理的變量變換方法，使變換后的“新因變量”與“新自變量”之間呈現(xiàn)出較好的相關(guān)關(guān)系，以便盡可能地反映圖1中大多數(shù)散點的變化趨勢(即具有較高的擬合優(yōu)度)。由此可知，“變量變換”是成功實現(xiàn)曲線擬合的有效途徑。

2 一個取自SAS幫助的數(shù)據(jù)集

2.1 數(shù)據(jù)集的名稱和數(shù)據(jù)結(jié)構(gòu)

在SAS幫助“數(shù)據(jù)庫”或“文件夾”中，有一個名為“sashelp.enso”的SAS數(shù)據(jù)集，其數(shù)據(jù)含義與結(jié)構(gòu)如表1所示[2]。

表1 澳大利亞港口城市與東部島嶼之間大氣壓值的變化數(shù)據(jù)

因篇幅限制，表1中僅列出了該數(shù)據(jù)集的前5行和最后5行。該數(shù)據(jù)集中含有3個變量和168個觀測(即N=168)。以下SAS程序可以顯示出完整的SAS數(shù)據(jù)集：

proc print data=sashelp.enso noobs;

run;

2.2 用散布圖呈現(xiàn)完整資料的變化趨勢

以下SAS程序可以呈現(xiàn)圖1中的大氣壓值(pressure)隨年份值(year)變化而變化的趨勢：

proc transreg data=sashelp.enso;

model identity(pressure)=identity(year);

run;

以上程序輸出結(jié)果為圖1。由圖1可知，除少數(shù)點外，絕大多數(shù)點幾乎是隨機地分布在一個長方形區(qū)域內(nèi)，近似呈“均勻分布”。此長方形與X軸不完全平行，略呈一個微小的傾斜角。圖1中的實線是基于最小二乘原理擬合出來的一條直線，按常規(guī)的統(tǒng)計學(xué)理念，這種表現(xiàn)的散布圖提示分析者，對此資料擬合直線回歸模型是無可非議的。然而，其R2僅為0.0071，說明用年份值去預(yù)測大氣壓值是非常不準(zhǔn)確的。

2.3 統(tǒng)計分析的任務(wù)

針對圖1中散點的分布或變化趨勢，可否找到一個較為合適的統(tǒng)計模型，使R2>0.6，即用年份值(year)預(yù)測大氣壓值(pressure)的預(yù)測結(jié)果具有一定程度的準(zhǔn)確性。

3 基于TRANSREG過程中各種樣條變換后建模[2-3]

3.1 基于B-樣條變換(spline)后建模

3.1.1 基本概念與做法

在SAS中，實現(xiàn)B-樣條變換的關(guān)鍵詞為“spline(自變量名)”。在運用SAS的TRANSREG過程時，可以對自變量year進行“B-樣條變換”，變換后的結(jié)果記為Tyear。再構(gòu)建因變量pressure關(guān)于新自變量Tyear的回歸模型。所謂“B-樣條變換”，實際上就是擬合因變量關(guān)于自變量的多項式曲線回歸模型，一次就是直線回歸模型、二次就是拋物線回歸模型、三次就是三次多項式回歸模型，以此類推。通常，擬合三次多項式回歸模型。

問題在于：是在自變量的整個取值區(qū)間內(nèi)擬合一個樣條函數(shù)曲線模型，還是將區(qū)間劃分成多個子區(qū)間，分別在各子區(qū)間上各擬合一個樣條函數(shù)曲線模型，即“分段擬合”。若在“nkonts=”之后寫一個“k(k≥0)”，就是將自變量的整個區(qū)間劃分成“k+1”個子區(qū)間?！皀konts”代表“節(jié)點數(shù)”或“斷點數(shù)”，“nkonts=k”代表節(jié)點數(shù)為“k”。顯然，隨著k值增加，分段數(shù)目也在增加，在各個很短的子區(qū)間上的“多項式曲線”能更好地擬合該子區(qū)間上的散點，因而擬合的效果就會提升，直到曲線回歸模型完全擬合給定的實際資料。下面給出“nkonts=0”“nkonts=5”“nkonts=10”等情形下的擬合結(jié)果。

3.1.2 SAS輸出結(jié)果及解釋

基于B-樣條變換且節(jié)點數(shù)k取不同數(shù)值時對應(yīng)的擬合效果(以R2來度量)。見表2。

表2 基于B-樣條變換且節(jié)點數(shù)k取不同數(shù)值時對應(yīng)的R2值

由表2可知，除了k=65和k=75兩行外，其他各行的R2都隨著k值增大而增大。其根本原因在前面已述及，此處不再贅述。按R2>0.6且回歸模型盡可能精簡的要求，k=35即可。若進一步嘗試，發(fā)現(xiàn)k=31時，R2=0.6751，此時擬合的圖形見圖2。

圖2 k=31時B-樣條變換對資料的擬合效果圖示

3.1.3 “nkonts=k”所對應(yīng)的SAS程序

注意：“k”必須取一個具體數(shù)值，下面給出k=31時對應(yīng)的SAS程序。

proc transreg data=sashelp.enso outtest=aaa ss2 coefficients test;

model identity(pressure)=spline(year/nknots=31);

output out=bbb predicted residuals;

run;

3.2 基于B-樣條基函數(shù)變換(bspline)后建模

3.2.1 基本概念與做法

在SAS中，實現(xiàn)B-樣條基函數(shù)變換的關(guān)鍵詞為“bspline(自變量名)”。與前面的“B-樣條變換”一樣，也可以通過給“nkonts=k”中的“k”賦一個具體值，來獲得擬合效果不同的擬合結(jié)果。而且，兩者的計算結(jié)果完全相同。它們之間的區(qū)別僅僅在于：使用“bspline(自變量名)”時，可以輸出“B-樣條基函數(shù)”的計算結(jié)果，而使用“spline(自變量名)”時，無法輸出“B-樣條基函數(shù)”的計算結(jié)果。

3.2.2 SAS輸出結(jié)果及解釋

基于B-樣條基函數(shù)變換且節(jié)點數(shù)k取不同數(shù)值時對應(yīng)的擬合效果(以R2來度量)，與表2完全相同，此處從略。下面分別將k=0、k=1和k=2時對應(yīng)的“B-樣條基函數(shù)”的計算結(jié)果呈現(xiàn)出來。

第一種情形：k=0時對應(yīng)的“B-樣條基函數(shù)”的計算結(jié)果見圖3。

圖3 k=0時對應(yīng)的“B-樣條基函數(shù)”的計算結(jié)果

圖3中的結(jié)果表明：在自變量year的整個取值區(qū)間內(nèi)，擬合了一個“三次多項式曲線回歸模型”，其表達式見式(1)。

(1)

第二種情形：k=1時對應(yīng)的“B-樣條基函數(shù)”的計算結(jié)果見圖4。

圖4 k=1時對應(yīng)的“B-樣條基函數(shù)”的計算結(jié)果

圖4中的結(jié)果表明：在前面式(1)基礎(chǔ)上，又增加了一項“t4”及其系數(shù)估計值。但前面各項的系數(shù)估計值也作了相應(yīng)的調(diào)整。

第三種情形：k=2時對應(yīng)的“B-樣條基函數(shù)”的計算結(jié)果見圖5。

圖5 k=2時對應(yīng)的“B-樣條基函數(shù)”的計算結(jié)果

圖5中的結(jié)果表明，在前面的基礎(chǔ)上，又增加了一項“t5”及其系數(shù)估計值。但前面各項的系數(shù)估計值也作了相應(yīng)的調(diào)整。

由此可知，自變量year經(jīng)過變量變換后的變量數(shù)目隨著k值的增加而增加，其具體個數(shù)為(4+k)。即k=0時，有4個；k=1時，有5個；而k=2時，有6個。值得注意的是，變換后的變量的下角標(biāo)從“0”開始。

3.2.3 “nkonts=k”所對應(yīng)的SAS程序

注意：“k”必須取一個具體數(shù)值，下面給出k=5時對應(yīng)的SAS程序。

proc transreg data=sashelp.enso outtest=aaa ss2 coefficients test;

model identity(pressure)=bspline(year/nknots=5);

output out=bbb predicted residuals;

run;

3.3 基于單調(diào)B-樣條變換后建模

3.3.1 基本概念與做法

在SAS中，實現(xiàn)單調(diào)B-樣條變換的關(guān)鍵詞為“mspline(自變量名)”。與前面的“B-樣條變換”一樣，也可以通過給“nkonts=k”中的“k”賦一個具體值，來獲得擬合效果不同的擬合結(jié)果。

3.3.2 SAS輸出結(jié)果及解釋

基于單調(diào)B-樣條變換且節(jié)點數(shù)k取不同數(shù)值時對應(yīng)的擬合效果(以R2來度量)。見表3。

表3 基于單調(diào)B-樣條變換且節(jié)點數(shù)k取不同數(shù)值時對應(yīng)的R2值

表3與表2相比，隨著回歸模型越來越復(fù)雜，表3中的R2值增加十分緩慢，且R2最大值均未超過0.05。

3.3.3 “nkonts=k”所對應(yīng)的SAS程序

注意：“k”必須取一個具體數(shù)值，下面給出k=160時對應(yīng)的SAS程序。

proc transreg data=sashelp.enso outtest=aaa ss2 coefficients test;

model identity(pressure)=mspline(year/nknots=160);

output out=bbb predicted residuals;

run;

3.4 基于非迭代懲罰B-樣條變換(pbspline)后建模

3.4.1 基本概念與做法

在SAS中，實現(xiàn)非迭代懲罰B-樣條變換的關(guān)鍵詞為“pbspline(自變量名)”。與前面的“B-樣條變換”一樣，也可以通過給“nkonts=k”中的“k”賦一個具體值，來獲得擬合效果不同的擬合結(jié)果。

3.4.2 SAS輸出結(jié)果及解釋

基于非迭代懲罰B-樣條變換且節(jié)點數(shù)k取不同數(shù)值時對應(yīng)的擬合效果(以R2來度量)。見表4。

表4 基于非迭代懲罰B-樣條變換且節(jié)點數(shù)k取不同數(shù)值時對應(yīng)的R2值

觀察表4，很難總結(jié)出其中規(guī)律，當(dāng)k=0、k=65、k≥100時，都能獲得較大的R2值。事實上，在使用pbspline變換時，涉及到一個關(guān)鍵的“光滑參數(shù)lambda”。當(dāng)此參數(shù)取不同數(shù)值時，回歸模型對資料的擬合效果是不同的。有多種評價指標(biāo)用來衡量回歸模型對資料的擬合優(yōu)度，其中，SBC統(tǒng)計量是經(jīng)常被選用的。SBC取值最小時，對應(yīng)的“l(fā)ambda值”是最合適的光滑參數(shù)。此時，對應(yīng)的R2將取極大值。

3.4.3 尋找使SBC取得極小值時所對應(yīng)的SAS程序

運行下面的一段SAS程序，可以找到使SBC=398.9為極小值，此時，lambda=1.14。

proc transreg data=sashelp.enso;

model identity(pressure) = pbspline(year/sbc lambda=.1 to 100.1 by 0.01);

run;

運行下面的一段SAS程序，可以找到使SBC=435.7為極小值，此時，lambda=1801.1。

proc transreg data=sashelp.enso;

model identity(pressure) = pbspline(year/sbc lambda=100.1 to 10000.1 by 0.01);

run;

由于在lambda的整個取值范圍內(nèi)，SBC取最小值時，回歸模型對資料的擬合優(yōu)度最好，故應(yīng)取lambda=1.14。使用下面的SAS程序，可以快速獲得最大的R2=0.6999。

proc transreg data=sashelp.enso outtest=aaa ss2 coefficients test;

model identity(pressure)=pbspline(year/sbc lambda=1.14);

output out=bbb predicted residuals;

run;

3.5 基于迭代光滑樣條變換(sspline)后建模

3.5.1 基本概念與做法

在SAS中，實現(xiàn)迭代光滑樣條變換的關(guān)鍵詞為“sspline(自變量名)”。與前面的“B-樣條變換”不一樣，它不能通過給“nkonts=k”中的“k”賦一個具體值，來獲得擬合效果不同的擬合結(jié)果。但它引入了一個“光滑參數(shù)smooth (簡寫為sm)”，通過調(diào)整sm的數(shù)值，可以獲得擬合效果不同的擬合結(jié)果。sm的取值從0開始，其取值區(qū)間為[0，100]，即可取此區(qū)間內(nèi)任何一個實數(shù)，也包括區(qū)間的兩個端點。數(shù)值越小，表明光滑程度越差，當(dāng)sm=0時，就是將所有相鄰的兩點用折線連起來，此時，為完全擬合(也稱為“過擬合”)；當(dāng)sm=100時，就是用一條直線來擬合該資料，擬合的效果最差。因此，可以依據(jù)R2的大小，嘗試尋找合適的sm數(shù)值。

3.5.2 SAS輸出結(jié)果及解釋

經(jīng)嘗試，當(dāng)sm=22.2時，R2=0.6039；當(dāng)sm=22.3時，R2=0.6008；當(dāng)sm=22.4時，R2=0.5976；當(dāng)sm=22.5時，R2=0.5944。通過進一步嘗試，發(fā)現(xiàn)當(dāng)sm=22.32時，R2=0.6001，剛好滿足R2>0.6的基本要求。

3.5.3尋找使R2略大于0.6時的sm數(shù)值所對應(yīng)的SAS程序

proc transreg data=sashelp.enso outtest=aaa ss2 coefficients test;

model identity(pressure)=sspline(year/sm=22.32);

output out=bbb predicted residuals;

run;

3.6 基于非迭代光滑樣條變換(smooth)后建模

3.6.1 基本概念與做法

在SAS中，實現(xiàn)非迭代光滑樣條變換的關(guān)鍵詞為“smooth(自變量名)”。與前面的“迭代光滑樣條變換”一樣，引入了一個“光滑參數(shù)smooth(簡寫成sm=)”，通過調(diào)整sm的數(shù)值，可以獲得擬合效果不同的擬合結(jié)果。

3.6.2 SAS輸出結(jié)果及解釋

經(jīng)嘗試，當(dāng)sm=20.8時，R2=0.6048；當(dāng)sm=20.9時，R2=0.6018；當(dāng)sm=21.0時，R2=0.5989；當(dāng)sm=21.1時，R2=0.5959。通過進一步嘗試，發(fā)現(xiàn)當(dāng)sm=20.96時，R2=0.6001，剛好滿足R2>0.6的基本要求。

3.6.3尋找使R2略大于0.6時的sm數(shù)值所對應(yīng)的SAS程序

proc transreg data=sashelp.enso outtest=aaa ss2 coefficients test;

model identity(pressure)=smooth(year/sm=20.96);

output out=bbb predicted residuals;

run;

4 討論與小結(jié)

本文圖1中的全部散點幾乎呈“均勻分布”狀態(tài)，常規(guī)的直線回歸模型沒有任何使用價值。然而，在SAS/STAT模塊中的“TRANSREG過程”收錄了很多種類的變量變換方法。本文采用了其中的一大類方法，即“樣條變換方法”。從思路和算法上又給出了六個彼此稍有區(qū)別的具體方法，它們各具優(yōu)點和缺點。根據(jù)分析者的目的不同，可以考慮選用不同的方法來擬合資料。

本文介紹的六種樣條變換方法主要作用是曲線擬合，在自變量的取值范圍內(nèi)找到合適的曲線回歸模型，即盡可能高的擬合優(yōu)度(如本文期望R2>0.6)且盡可能精簡的回歸模型(可以回歸模型誤差項的自由度“DF誤差”來反映，其數(shù)值越大越好，它標(biāo)志著回歸模型中被估計的參數(shù)數(shù)目少，即回歸模型較為精簡)。

B-樣條變換和B-樣條基函數(shù)變換：節(jié)點數(shù)為31時，R2=0.6751，DF誤差=133。

單調(diào)B-樣條變換：擬合效果極差。

非迭代懲罰B-樣條變換：Lambda=1.14時，R2=0.6999,DF誤差=131.28。

迭代光滑樣條變換：Sm=22.32時，R2=0.6001，DF誤差=138.34。

非迭代光滑樣條變換：Sm=20.96時，R2=0.6001，DF誤差=137.09。

由于以上方法所依賴的“調(diào)節(jié)參數(shù)”不同，因此，結(jié)果尚缺乏可比性。但可將R2值定為0.6999，計算出除“單調(diào)B-樣條變換”之外的其他幾種情況下的結(jié)果：

B-樣條變換和B-樣條基函數(shù)變換：節(jié)點數(shù)為34時，R2=0.6974，DF誤差=130；節(jié)點數(shù)為35時，R2=0.7171，DF誤差=129。

非迭代懲罰B-樣條變換：Lambda=1.14時，R2=0.6999,DF誤差=131.28。

迭代光滑樣條變換：Sm=18.4時，R2=0.6998，DF誤差=131.34。

非迭代光滑樣條變換：Sm=16.88時，R2=0.6999，DF誤差=129.18。

綜上所述，“非迭代懲罰B-樣條變換”與“迭代光滑樣條變換”的擬合效果基本相同，是以上六種方法中最好的曲線回歸建模方法。