變量變換回歸分析(II)擬合近似呈均勻分布資料的方法
——

胡良平

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

1 概 述

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

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

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

1.2 變量變換的必要性

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

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

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

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

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

因篇幅限制，表1中僅列出了該數(shù)據(jù)集的前5行和最后5行。該數(shù)據(jù)集中含有3個(gè)變量和168個(gè)觀測(即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ù)點(diǎn)外，絕大多數(shù)點(diǎn)幾乎是隨機(jī)地分布在一個(gè)長方形區(qū)域內(nèi)，近似呈“均勻分布”。此長方形與X軸不完全平行，略呈一個(gè)微小的傾斜角。圖1中的實(shí)線是基于最小二乘原理擬合出來的一條直線，按常規(guī)的統(tǒng)計(jì)學(xué)理念，這種表現(xiàn)的散布圖提示分析者，對此資料擬合直線回歸模型是無可非議的。然而，其R2僅為0.0071，說明用年份值去預(yù)測大氣壓值是非常不準(zhǔn)確的。

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

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

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

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

3.1.1 基本概念與做法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

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

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

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

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

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

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

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

注意：“k”必須取一個(gè)具體數(shù)值，下面給出k=5時(shí)對應(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中，實(shí)現(xiàn)單調(diào)B-樣條變換的關(guān)鍵詞為“mspline(自變量名)”。與前面的“B-樣條變換”一樣，也可以通過給“nkonts=k”中的“k”賦一個(gè)具體值，來獲得擬合效果不同的擬合結(jié)果。

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

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

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

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

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

注意：“k”必須取一個(gè)具體數(shù)值，下面給出k=160時(shí)對應(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中，實(shí)現(xiàn)非迭代懲罰B-樣條變換的關(guān)鍵詞為“pbspline(自變量名)”。與前面的“B-樣條變換”一樣，也可以通過給“nkonts=k”中的“k”賦一個(gè)具體值，來獲得擬合效果不同的擬合結(jié)果。

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

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

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

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

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

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

proc transreg data=sashelp.enso;

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

run;

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

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

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

3.5.3尋找使R2略大于0.6時(shí)的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中，實(shí)現(xiàn)非迭代光滑樣條變換的關(guān)鍵詞為“smooth(自變量名)”。與前面的“迭代光滑樣條變換”一樣，引入了一個(gè)“光滑參數(shù)smooth(簡寫成sm=)”，通過調(diào)整sm的數(shù)值，可以獲得擬合效果不同的擬合結(jié)果。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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