胡良平

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

1 概 述[1]

1.1 局部回歸模型

局部回歸模型見式(1)：

yi=g(xi)+εii=1,2,…,n

(1)

在式(1)中，yi為第i次觀測到的因變量的取值；g(xi)是xi的回歸函數(shù)；xi可以是一個(gè)自變量，也可以是由多個(gè)自變量組成的向量；εi是一個(gè)隨機(jī)誤差。

1.2 局部模型回歸分析應(yīng)用的場合

一般來說，在因變量服從正態(tài)分布或?qū)ΨQ分布時(shí)，欲研究因變量隨自變量變化而變化的依賴關(guān)系時(shí)，可以嘗試采用很多種方法來創(chuàng)建回歸模型，包括采用“局部回歸模型”。最適合運(yùn)用此模型的場合如下：在自變量的全部取值范圍內(nèi)，存在多個(gè)“小區(qū)域”，在這些“小區(qū)域”內(nèi)，觀測點(diǎn)的密度較高，似乎呈現(xiàn)出“聚集性”；而且，它們或呈“二次多項(xiàng)式曲線形狀”或呈“三次多項(xiàng)式曲線形狀”分布。見圖1。

圖1 黑色素瘤發(fā)病率隨時(shí)間推移的變化趨勢

1.3 局部模型回歸分析的計(jì)算原理

1.3.1 計(jì)算原理

所謂局部模型，實(shí)際上就是在每個(gè)“小區(qū)域或小鄰域”上構(gòu)建自變量的一個(gè)線性或二次曲線模型、甚至三次曲線模型。問題在于如何選取一系列的“小鄰域”。一個(gè)最直觀的想法是：將全部數(shù)據(jù)觀察點(diǎn)按自變量由小到大的順序排列，先確定由多少個(gè)相鄰的觀察點(diǎn)決定一個(gè)“小鄰域”，比如，設(shè)觀察點(diǎn)數(shù)目為k(k≥3)，當(dāng)k取一個(gè)確定數(shù)值后，就很容易將全部觀察點(diǎn)劃分成m個(gè)“小鄰域”。于是，在每個(gè)“小鄰域”上創(chuàng)建一個(gè)“局部模型”，計(jì)算出各“小鄰域”上因變量的殘差平方和，再求出所有“小鄰域”上殘差平方法和之和，就可獲得總殘差平方和。接下去，就可以改變k值，假定令k=3到k=n(即全部觀察點(diǎn))共有j種情況，由前面的計(jì)算就可獲得某種情況下的“總殘差平方和”最小，于是，就認(rèn)為按這種情況對(duì)應(yīng)的“k值”來形成“小鄰域”是最合適的。

事實(shí)上，在SAS的LOESS過程中，評(píng)價(jià)擬合效果所選用的統(tǒng)計(jì)量為校正的赤池信息準(zhǔn)則(AICC )(其取值越小越好，具體計(jì)算公式詳見后文)，它所對(duì)應(yīng)的k值被轉(zhuǎn)換成“光滑參數(shù)s”，s=k/n(其中k需要事先依據(jù)某種方法或理由初步估計(jì)出來，n為樣本含量或全部觀察點(diǎn)數(shù)目)。在每個(gè)“小鄰域”上建模時(shí)，采用“加權(quán)最小平方法”[2]。

1.3.2 常用的擬合效果評(píng)價(jià)指標(biāo)

(1)赤池信息準(zhǔn)則(The Akaike information criterion，AIC)：AIC是模型對(duì)資料擬合優(yōu)度的一種度量，也體現(xiàn)了現(xiàn)在所使用的模型相對(duì)于最簡約模型之間的一種平衡。其定義如下：

AIC=-2LL+2p

上式中，p為模型中被估計(jì)參數(shù)的個(gè)數(shù)，LL是用于估計(jì)參數(shù)數(shù)值的似然函數(shù)的對(duì)數(shù)。

(2)AICC：

上式中，n為總樣本含量，其他變量含義同上。

(3)貝葉斯信息準(zhǔn)則(Bayesian Information Criterion，BIC)與AIC和AICC是類似的度量，其定義如下：

BIC=-2LL+plog(n)

上式中，各變量的含義同上，此處不再贅述。

2 基于局部模型回歸分析解決實(shí)際問題[1]

2.1 問題與數(shù)據(jù)結(jié)構(gòu)

【例1】下面是一個(gè)關(guān)于黑色素瘤發(fā)病率的資料。資料來自美國康涅狄格州腫瘤注冊(cè)部門，時(shí)間從1936年-1972年共37年，基于年齡校正的各年黑色素瘤的發(fā)病率(1/10萬)的前8年數(shù)據(jù)見表1，其他數(shù)據(jù)詳見后面的SAS程序：

表1 基于年齡校正的1936年-1943年黑色素瘤發(fā)病率

【對(duì)數(shù)據(jù)結(jié)構(gòu)的分析】嚴(yán)格地說，這是一個(gè)“時(shí)間序列”數(shù)據(jù)，即發(fā)病率隨著時(shí)間的推移而動(dòng)態(tài)變化。為簡便起見，暫且將該數(shù)據(jù)視為一個(gè)計(jì)量因變量y(發(fā)病率)隨另一個(gè)計(jì)量自變量x(年份)變化的依賴關(guān)系問題。

【統(tǒng)計(jì)分析方法的選擇】研究y與x之間依賴關(guān)系的最簡單方法是進(jìn)行直線回歸分析；若兩變量之間呈曲線變化趨勢，就可選擇某種曲線方程進(jìn)行曲線回歸分析。

2.2 基于常規(guī)方法構(gòu)建簡單線性回歸模型[3]

2.2.1 創(chuàng)建SAS數(shù)據(jù)集

創(chuàng)建一個(gè)名為“melanoma”的臨時(shí)SAS數(shù)據(jù)集的SAS數(shù)據(jù)步程序如下：

data Melanoma;

input Year Incidences @@;

format Year d4.0;

datalines;

19360.919370.819380.819391.319401.419411.219421.719431.819441.619451.519461.519472.019482.519492.719502.919512.519523.119532.419542.219552.919562.519572.619583.219593.819604.219613.919623.719633.319643.719653.919664.119673.819684.719694.419704.819714.819724.8

;

run;

2.2.2 繪制散布圖，直觀展示兩變量之間的變化趨勢

利用下面的SAS過程步程序，可以繪制反映兩變量變化趨勢：

proc sgplot data=Melanoma;

scatter y=Incidences x=Year;

run;

【SAS輸出結(jié)果】

第1部分輸出結(jié)果為“圖1”，已經(jīng)在前面呈現(xiàn)，此處從略。

由圖1可看出：散點(diǎn)呈上升的變化趨勢。但仔細(xì)觀察散點(diǎn)，發(fā)現(xiàn)在多個(gè)局部區(qū)域內(nèi)散點(diǎn)表現(xiàn)為“聚集性”，并且呈“矩形”或“三角形”等形狀。

下面嘗試采用簡單直線回歸模型擬合該資料：

ods graphics on;

proc reg data=Melanoma;

model Incidences=Year;

run;

【SAS主要輸出結(jié)果】

圖2 采用直線回歸模型描述黑色素瘤發(fā)病率隨時(shí)間推移的變化趨勢

擬合的統(tǒng)計(jì)量：均方根誤差=0.33641、R2=0.9283、調(diào)整R2=0.9263，從這些擬合統(tǒng)計(jì)量的數(shù)值來看，似乎用簡單直線回歸模型擬合此資料效果相當(dāng)令人滿意。但從圖2可看出：在多個(gè)局部區(qū)域上，直線不能很好地給出預(yù)測結(jié)果。

2.3 基于局部模型構(gòu)建非線性回歸模型[1]

基于局部模型構(gòu)建非線性回歸模型的SAS程序如下：

proc loess data=Melanoma;

model Incidences=Year;

run;

【SAS程序說明】以上SAS程序調(diào)用LOESS過程擬合局部模型。

【SAS輸出結(jié)果及其解釋】

由圖3可看出：局部模型對(duì)此資料的擬合效果非常好，既沒有“過擬合”，也沒有“欠擬合”。

如何才能做到既不“過擬合”又不“欠擬合”？關(guān)鍵是要選取合適的“光滑參數(shù)”，它已顯示在圖3的左上角，即“Smooth=0.257”。用此數(shù)值乘以總樣本含量37等于9.5，說明程序按橫坐標(biāo)軸的順序，將每相鄰9或10個(gè)觀測點(diǎn)所在的區(qū)域視為一個(gè)“局部區(qū)域”，在該區(qū)域上進(jìn)行多項(xiàng)式擬合。

圖3 采用局部模型擬合的結(jié)果

如何獲得最佳“光滑參數(shù)”的數(shù)值？在SAS的LOESS過程中，先給定一系列的“光滑參數(shù)”值進(jìn)行擬合，對(duì)于每個(gè)給定的“光滑參數(shù)”值，就能計(jì)算出若干個(gè)反映擬合效果或優(yōu)度的統(tǒng)計(jì)量，其中，以AICC統(tǒng)計(jì)量取得最小值時(shí)對(duì)應(yīng)的“光滑參數(shù)”為最佳。

利用如下SAS程序可以同時(shí)獲得4個(gè)“光滑參數(shù)”對(duì)應(yīng)的擬合結(jié)果，

proc loess data=Melanoma plots=ResidualsBySmooth(smooth);

model Incidences=Year/smooth=0.1 0.25 0.4 0.6;

run;

【SAS主要輸出結(jié)果】

圖4 基于4個(gè)光滑參數(shù)進(jìn)行局部模型擬合得到的擬合結(jié)果

在圖4中有4幅小圖，從上往下、從左往右的“光滑參數(shù)”依次為0.1、0.25、0.4和0.6對(duì)應(yīng)的擬合結(jié)果。不難看出：“Smooth=0.1”屬于“過擬合”，而“Smooth=0.4”和“Smooth=0.6”屬于“欠擬合”，只有“Smooth=0.25”，屬于“正常擬合”，因?yàn)樗呀?jīng)是最佳“光滑參數(shù)”0.257的近似值。

圖5 基于4個(gè)光滑參數(shù)進(jìn)行局部模型擬合得到的殘差圖

圖5中的4幅小圖分別與圖4中4幅小圖一一對(duì)應(yīng)，只不過圖5反映的是殘差。當(dāng)“Smooth=0.1”時(shí)，幾乎所有觀察點(diǎn)上的殘差都為0，這就是“過擬合”；當(dāng)“Smooth=0.25”時(shí)，殘差圖上散點(diǎn)在各處波動(dòng)接近且沒有明顯的變化趨勢，屬于“正常擬合”；而圖5中下面的2幅小圖都呈現(xiàn)出殘差散點(diǎn)具有一定的變化規(guī)律，屬于“欠擬合”。

為了避免盲目性，可以采用下面的SAS程序自動(dòng)尋找到最佳的“光滑參數(shù)”的數(shù)值：

proc loess data=Melanoma;

model Incidences=Year / details(ModelSummary OutputStatistics);

run;

【SAS主要輸出結(jié)果】

Model SummarySmoothingParameterLocalPointsResidual SSGCVAICC0.41892153.422290.00339-0.962520.68919254.058380.00359-0.934590.31081112.510540.00279-1.120340.2027071.585130.00239-1.122210.1756861.568960.00241-1.097060.28378102.504870.00282-1.104020.2027071.585130.00239-1.122210.2567692.031050.00252-1.172770.2297382.029650.00256-1.151450.2567692.031050.00252-1.17277

以上是程序自動(dòng)尋找最佳“光滑參數(shù)”的動(dòng)態(tài)過程，僅當(dāng)局部觀測點(diǎn)為9個(gè)時(shí)，AICC統(tǒng)計(jì)量能取到最小值-1.17277，此時(shí)，對(duì)應(yīng)的“光滑參數(shù)”為0.25676。

Fit SummaryFit Methodkd TreeBlendingLinearNumber of Observations37Number of Fitting Points37kd Tree Bucket Size1Degree of Local Polynomials1Smoothing Parameter0.25676Points in Local Neighborhood9Residual Sum of Squares2.03105Trace[L]8.62243GCV0.00252AICC-1.17277

以上是模型擬合效果的總結(jié)。

利用下面的SAS程序，可以得到擬合曲線的置信帶：

proc loess data=Melanoma;

model Incidences=Year/clm alpha=0.05;

run;

【SAS主要輸出結(jié)果】

Fit SummaryFit Methodkd TreeBlendingLinearNumber of Observations37Number of Fitting Points37kd Tree Bucket Size1Degree of Local Polynomials1Smoothing Parameter0.25676Points in Local Neighborhood9Residual Sum of Squares2.03105Trace[L]8.62243GCV0.00252AICC-1.17277AICC1-42.03789Delta127.06596Delta226.76564Equivalent Number of Parameters7.31083Lookup Degrees of Freedom27.36964Residual Standard Error0.27394

以上是模型擬合效果的總結(jié)，與前面給出的結(jié)果基本相同。

圖6 基于光滑參數(shù)為0.257時(shí)得到的局部多項(xiàng)式擬合結(jié)果及95%置信帶

2.4 小結(jié)

從上面的介紹可知：局部模型的關(guān)鍵在于選取“光滑參數(shù)”的具體取值。此值的真實(shí)含義是以每相鄰的多少個(gè)觀察點(diǎn)為一個(gè)“小區(qū)域”，在每個(gè)這樣的“小區(qū)域”上擬合一個(gè)“多項(xiàng)式”。當(dāng)“Smooth=0.1”(相當(dāng)于樣本含量的1/10的觀察點(diǎn))時(shí)，得到了“過擬合”的結(jié)果。就本例而言，37/10=3.7≈4，若采用4次多項(xiàng)式，則多項(xiàng)式曲線就會(huì)通過每個(gè)觀察點(diǎn)；當(dāng)“Smooth=0.6”(相當(dāng)于樣本含量的6/10的觀察點(diǎn))時(shí)，得到了“欠擬合”的結(jié)果。就本例而言，6×(37/10)≈22，若采用4次多項(xiàng)式，則多項(xiàng)式曲線就很難通過大多數(shù)觀察點(diǎn)。

當(dāng)采用簡單直線回歸模型時(shí)，就相當(dāng)于取“Smooth=1.0”，也就把全部觀察點(diǎn)所在的范圍視為一個(gè)“小區(qū)域”，采用一個(gè)“一次多項(xiàng)式”去擬合資料，這對(duì)于具有類似圖1中散點(diǎn)所表現(xiàn)的狀態(tài)是沒有任何幫助的。

由此可知：局部模型最適合用于如下的資料：全部觀察點(diǎn)呈現(xiàn)線性遞增或下降趨勢，而在多個(gè)“小區(qū)域”上表現(xiàn)為“二次曲線”或“三次曲線”或“四次曲線”的形狀。建模的目的只是為了形象化地?cái)M合數(shù)據(jù)并對(duì)未知因變量的取值進(jìn)行預(yù)測，而不需要呈現(xiàn)回歸模型的具體表達(dá)式(因此法不便給出具體的回歸模型)。