胡良平

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

1 基本概念

1.1 變量變換

在回歸分析中，常需要對(duì)變量進(jìn)行變換，以使資料滿足待擬合回歸模型的要求。最常見(jiàn)的情形是對(duì)定性自變量作啞變量變換，對(duì)定量變量作對(duì)數(shù)變換、指數(shù)變換、平方根變換或平方變換等。根據(jù)變量變換所涉及的原理和計(jì)算方法的不同，在SAS的TRANSREG過(guò)程中，變量變換方法分為五大類，每一類中又可以細(xì)分成許多種小類[1]。因篇幅所限，此處暫不贅述。

1.2 樣條變換

“樣條變換”是變量變換中頗具特色的一類變換方法，其基本思路是擬合“分段多項(xiàng)式曲線”。若根據(jù)具體算法，樣條變換可以進(jìn)一步劃分為“B-樣條變換”“B-樣條基函數(shù)變換”“單調(diào)B-樣條變換”“無(wú)迭代懲罰B-樣條變換”“迭代光滑樣條變換”和“非迭代光滑樣條變換”等子類[1]。

1.3 節(jié)點(diǎn)

所謂“節(jié)點(diǎn)”，就是“間斷點(diǎn)”或“不連續(xù)點(diǎn)”。通常有兩種“節(jié)點(diǎn)”：其一，實(shí)際資料中客觀存在的“間斷點(diǎn)”，例如本文的圖1中就出現(xiàn)了3個(gè)“節(jié)點(diǎn)”；其二，統(tǒng)計(jì)分析者人為設(shè)定的“間斷點(diǎn)”，例如根據(jù)某實(shí)際資料的散布圖中各散點(diǎn)的分布情況，似乎可以分別采用“對(duì)數(shù)曲線”“指數(shù)曲線”和“冪函數(shù)曲線”去描述該實(shí)際資料的變化趨勢(shì)，于是，分析者就需要在自變量的取值區(qū)間內(nèi)設(shè)定兩個(gè)“節(jié)點(diǎn)”，將整個(gè)曲線劃分成三段。對(duì)此資料的擬合優(yōu)度高低，取決于兩個(gè)“節(jié)點(diǎn)”的位置選得是否恰當(dāng)。本文將針對(duì)前述提及的兩種“節(jié)點(diǎn)”分別擬合曲線，以說(shuō)明在進(jìn)行曲線擬合時(shí)，如何發(fā)揮“節(jié)點(diǎn)”的作用。

2 一個(gè)人工生成的實(shí)例——由三個(gè)間斷點(diǎn)分割的四段曲線

2.1 產(chǎn)生實(shí)例數(shù)據(jù)集的SAS程序

利用下面的SAS程序，可以生成一個(gè)由三個(gè)間斷點(diǎn)分割的四段曲線。

title1'An Illustration of Splines and Knots';

* Create in y a discontinuous function of x.;

data a;

x=-0.000001;

do i=0 to 199;

if mod(i, 50) = 0 then do;

c=((x/2)-5)**2;

if i=150 then c=c+5;

y=c;

end;

x=x+0.1;

y=y-sin(x-c);

output;

end;

run;

2.2 數(shù)據(jù)集a的基本情況描述與呈現(xiàn)

數(shù)據(jù)集a中主要有兩個(gè)計(jì)量變量，即x與y，共有200個(gè)觀測(cè)點(diǎn)。x與y的簡(jiǎn)單統(tǒng)計(jì)量如下。

變量N均值標(biāo)準(zhǔn)差最小值最大值x20010.04999905.78791850.099999019.9999990y20012.04335399.2339544-7.370567029.1657155

其中，x呈均勻分布、y呈負(fù)偏態(tài)分布(圖形從略)。(x，y)的散布圖見(jiàn)圖1。

圖1 基于數(shù)據(jù)集a繪制出(x，y)的散布圖

由圖1可知，在變量x的整個(gè)變化范圍[0.0999990,19.9999990]中，共有4段曲線，它們被3個(gè)不連續(xù)點(diǎn)(x=5、10、15)分隔開(kāi)。從左至右，第一、三和四段曲線很像二次拋物線，第二段曲線很像三次拋物線。

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

試以x為自變量、y為因變量，建立y依賴x變化而變化的回歸模型。換言之，要對(duì)圖1所代表的數(shù)據(jù)集進(jìn)行曲線擬合，希望構(gòu)建一個(gè)回歸方程或回歸模型，以便用區(qū)間[0.0999990,19.9999990]中的任何一個(gè)x的取值代入，求出因變量y的估計(jì)值。要求是預(yù)測(cè)誤差應(yīng)盡可能小。

3 不考慮與考慮節(jié)點(diǎn)且使用多項(xiàng)式函數(shù)進(jìn)行曲線擬合

3.1 不考慮節(jié)點(diǎn)且僅使用多項(xiàng)式函數(shù)進(jìn)行曲線擬合

3.1.1 呈現(xiàn)各種情況下的擬合效果

所謂“不考慮節(jié)點(diǎn)”，就是將全部數(shù)據(jù)集視為一個(gè)“整體”，擬合一條光滑的曲線回歸方程或模型。通常，可以考慮從擬合2次多項(xiàng)式開(kāi)始，逐漸增加“次數(shù)”。以下依次擬合2～16次多項(xiàng)式，以“R2”反映其擬合優(yōu)度。擬合效果見(jiàn)表1。

表1 對(duì)數(shù)據(jù)集a擬合2～16次多項(xiàng)式所得的擬合效果

由表1可知，隨著多項(xiàng)式“次數(shù)”的增加，R2也不斷增大。但是，多項(xiàng)式的次數(shù)每增加1次，就相當(dāng)于增加了一個(gè)新變量。一個(gè)包含16次多項(xiàng)式的回歸模型中包含了17個(gè)未知參數(shù)(因?yàn)檫€包含1個(gè)常數(shù)項(xiàng))，顯然增加了模型的復(fù)雜程度，模型的實(shí)用性大大降低了。不僅如此，模型在“不連續(xù)點(diǎn)(x=5、10、15)”附近的擬合誤差很大。見(jiàn)圖2。

圖2 基于樣條變換且采用16次多項(xiàng)式擬合數(shù)據(jù)集a的效果

3.1.2多項(xiàng)式次數(shù)為16時(shí)對(duì)應(yīng)的SAS過(guò)程步程序

表1中最后一行(多項(xiàng)式的次數(shù)為16)所對(duì)應(yīng)的SAS程序如下：

proc transreg data=a;

model identity(y)=spline(x / degree=16);

run;quit;

依次將上述程序中的選項(xiàng)“degree=”賦值為2～15，可獲得表1中其他各行的計(jì)算結(jié)果。

3.2 由程序確定節(jié)點(diǎn)位置且使用樣條變換進(jìn)行曲線擬合

3.2.1 呈現(xiàn)各種情況下的擬合效果

分析者可以指定節(jié)點(diǎn)的個(gè)數(shù)，一旦個(gè)數(shù)確定，將由程序按照規(guī)定的方式計(jì)算出各節(jié)點(diǎn)所在的位置。于是，整個(gè)數(shù)據(jù)集就被分割為若干個(gè)子區(qū)間，程序?qū)⒃诿總€(gè)子區(qū)間上擬合樣條函數(shù)，這就是所謂的“分段擬合”，分段的數(shù)目為節(jié)點(diǎn)數(shù)加1。在每一段上，又可以擬合2、3、4等高次多項(xiàng)式。但通常的樣條函數(shù)為3次多項(xiàng)式，即分段擬合3次多項(xiàng)式。下面僅改變節(jié)點(diǎn)數(shù)，但在每段上都將擬合3次多項(xiàng)式。節(jié)點(diǎn)數(shù)從1到15，對(duì)應(yīng)的擬合效果見(jiàn)表2。

由表2可知，從總趨勢(shì)上來(lái)看，隨著節(jié)點(diǎn)數(shù)增加，R2也在增大，但當(dāng)節(jié)點(diǎn)數(shù)為5、8和11時(shí)，R2略有下降。由于節(jié)點(diǎn)數(shù)越多，整個(gè)區(qū)間被劃分成的子區(qū)間數(shù)目就越多，也即3次多項(xiàng)式的個(gè)數(shù)越多，回歸模型越復(fù)雜。下面給出當(dāng)R2取得最大值0.97508(節(jié)點(diǎn)數(shù)為14)時(shí)，對(duì)應(yīng)的曲線擬合結(jié)果。見(jiàn)圖3。

表2考慮不同節(jié)點(diǎn)數(shù)且僅擬合3次多項(xiàng)式對(duì)數(shù)據(jù)集a擬合的效果

編 號(hào)節(jié)點(diǎn)數(shù)R2編 號(hào)節(jié)點(diǎn)數(shù)R2110.53585990.95256220.5539110100.96674330.6146511110.95835440.7445012120.96231550.6882313130.97174660.8889014140.97508770.9360815150.96872880.92566

圖3 基于節(jié)點(diǎn)數(shù)為14且采用3次樣條擬合數(shù)據(jù)集a的效果

3.2.2 節(jié)點(diǎn)數(shù)為15時(shí)對(duì)應(yīng)的SAS過(guò)程步程序

表2中最后一行所對(duì)應(yīng)的SAS程序如下：

proc transreg data=a;

model identity(y) = spline(x/degree=3 nknots=15);

run;

依次將上述程序中的選項(xiàng)“nknots=”賦值為1～14，可獲得表2中其他各行的計(jì)算結(jié)果。

3.3 基于數(shù)據(jù)集的真實(shí)情況給定3個(gè)節(jié)點(diǎn)位置且使用3次樣條變換進(jìn)行曲線擬合

3.3.1每個(gè)節(jié)點(diǎn)出現(xiàn)k(k=1～4)次

在x=5、10、15三處真實(shí)的“不連續(xù)點(diǎn)”處，對(duì)函數(shù)求一階導(dǎo)數(shù)；若將三個(gè)不連續(xù)點(diǎn)各重復(fù)寫2次，就是對(duì)函數(shù)求二階導(dǎo)數(shù)；同理，若將三個(gè)不連續(xù)點(diǎn)各重復(fù)寫k(k=1～10)次，就是對(duì)函數(shù)求k階導(dǎo)數(shù)。求不同階導(dǎo)數(shù)后，擬合的效果可能會(huì)有所變化。經(jīng)嘗試，當(dāng)k≥4時(shí)，擬合的效果幾乎相同。見(jiàn)表3。

表3 將三個(gè)“不連續(xù)點(diǎn)”各重復(fù)k次且僅用3次多項(xiàng)式對(duì)數(shù)據(jù)集a擬合的效果

由表3可知，將三個(gè)“不連續(xù)點(diǎn)”各重復(fù)4次，就能獲得最好的擬合效果。見(jiàn)圖4。

圖4 將三個(gè)“不連續(xù)點(diǎn)”各重復(fù)4次且僅用3次多項(xiàng)式的擬合效果

3.3.2將三個(gè)“不連續(xù)點(diǎn)”各重復(fù)4次時(shí)對(duì)應(yīng)的SAS過(guò)程步程序

proc transreg data=a;

model identity(y) = spline(x / knots=

5 5 5 5

10 10 10 10

15 15 15 15);

run;quit;

4 討論與小結(jié)

常規(guī)的曲線擬合方法不適合含有間斷點(diǎn)的資料，其主要原因在于：常規(guī)的曲線類型基本上都在自變量的定義域范圍內(nèi)是連續(xù)的[2-4]。本文數(shù)據(jù)集a含有3個(gè)間斷點(diǎn)，故常規(guī)的曲線擬合方法不適合。

本文通過(guò)三種方法來(lái)擬合含有間斷點(diǎn)的資料。第一種：不考慮節(jié)點(diǎn)且僅考慮多項(xiàng)式的次數(shù)，由表1可知，欲使R2達(dá)到0.96以上，至少需要采用14次多項(xiàng)式；第二種：基于3次多項(xiàng)式且考慮節(jié)點(diǎn)的數(shù)目，由表2可知，欲使R2達(dá)到0.96以上，至少需要采用10個(gè)節(jié)點(diǎn)；第三種：基于3次多項(xiàng)式、瞄準(zhǔn)資料中真實(shí)的3個(gè)“間斷點(diǎn)”，并強(qiáng)調(diào)在“間斷點(diǎn)”處重復(fù)不同次數(shù)，欲使R2達(dá)到0.96以上，只需要在3個(gè)“間斷點(diǎn)”處重復(fù)4次即可。

欲達(dá)到相同或相近的擬合效果，第一種方法所對(duì)應(yīng)的回歸模型過(guò)于復(fù)雜且實(shí)用性很差；第二種方法所對(duì)應(yīng)的回歸模型有了很大的改進(jìn)，回歸模型的實(shí)用性有了較大提升；第三種方法所對(duì)應(yīng)的回歸模型最簡(jiǎn)潔，也最具有實(shí)用性。