朋文佳，朱玉，賈賢杰，羋靜

科學(xué)研究的一個(gè)重要目的是評估結(jié)局變量與解釋變量之間的數(shù)量依存關(guān)系，回歸分析是達(dá)到該研究目的重要統(tǒng)計(jì)學(xué)分析方法?；貧w分析存在著重要價(jià)值，其一描述結(jié)局變量與解釋變量之間的數(shù)量依存關(guān)系；其二通過解釋變量預(yù)測結(jié)局變量；其三比較眾多解釋變量對結(jié)局變量影響大小[1]?，F(xiàn)實(shí)世界中，結(jié)局變量與解釋變量的數(shù)量依存關(guān)系錯綜復(fù)雜，為了更好地反映數(shù)量依存關(guān)系，回歸分析理論在逐步完善，形成了參數(shù)回歸分析和非參數(shù)回歸分析[2]。在教學(xué)中發(fā)現(xiàn)學(xué)生對參數(shù)回歸和非參數(shù)回歸的理解和應(yīng)用存在困難，本文以科研案例為基礎(chǔ)，將參數(shù)回歸分析和非參數(shù)回歸分析的關(guān)鍵理論進(jìn)行總結(jié)和辨析，并用R語言實(shí)現(xiàn)分析。

1 案例數(shù)據(jù)與分析目的

1992-1996年期間，美國開展了一項(xiàng)“脊柱骨密度”科學(xué)研究[3]，本文從該研究中截取了部分女性數(shù)據(jù)，見表1。數(shù)據(jù)集包含423位調(diào)查對象，3個(gè)變量。其中，變量id為調(diào)查對象編號；變量y是調(diào)查對象的脊柱骨密度含量，為結(jié)局變量；變量x是調(diào)查對象的年齡，為解釋變量。本案例欲研究年齡對脊柱骨密度含量有無影響？如果存在影響，影響形式是什么？

表1 脊柱骨密度研究數(shù)據(jù)

2 參數(shù)回歸及R語言實(shí)現(xiàn)

2.1 參數(shù)回歸簡介 因?yàn)榧怪敲芏群渴沁B續(xù)型定量變量，則解釋變量年齡(x)與結(jié)局變量脊柱骨密度(y)之間的潛在趨勢用參數(shù)回歸形式描述為：

yi=f(xi) +εi,E(yi|xi)=f(xi)=β0+β1xi,E(εi)=0

其中函數(shù)f(xi) =β0+β1xi形式固定，β0和β1為待估計(jì)的未知參數(shù)，εi為隨機(jī)誤差，i=1,2,L,n為個(gè)體編號。因?yàn)閒(xi) 是固定的函數(shù)形式，函數(shù)中未知參數(shù)有明確的統(tǒng)計(jì)含義，該分析方法被稱為線性回歸分析，是最簡單的參數(shù)回歸。β0和β1可以利用最小二乘法估計(jì)，即殘差平和(SS)最小[4-5]：

線性回歸要求數(shù)據(jù)滿足“線性”、“獨(dú)立”、“正態(tài)分布”和“方差齊”四個(gè)條件。線性回歸中函數(shù)左邊是y均數(shù)恒等變換，右邊是x的線性表達(dá)式，所以在直角坐標(biāo)系中數(shù)據(jù)點(diǎn)(x，y)呈現(xiàn)出“線性”；“獨(dú)立”要求每個(gè)數(shù)據(jù)點(diǎn)相互獨(dú)立；“正態(tài)分布”是指任意x取值下，y分布是正態(tài)分布；“方差齊”是指任意x取值下，y的方差保持不變。線性回歸的要求限制了線性回歸的使用，為了滿足更多數(shù)據(jù)分析需求，將結(jié)局變量的分布擴(kuò)展到指數(shù)分布族，形成了廣義線性回歸[6]。進(jìn)行線性回歸分析可以分為三階段，其一繪制散點(diǎn)圖，其二最小二乘法估計(jì)參數(shù)和統(tǒng)計(jì)推斷，其三回歸診斷。

2.2 參數(shù)回歸R語言實(shí)現(xiàn) 繪制散點(diǎn)圖(見圖1)，可見隨著年齡增長，脊柱骨密度含量上升，年齡前段脊柱骨密度含量增長多，年齡后段脊柱骨密度含量增長少。R語言是非常流行的統(tǒng)計(jì)分析語言。在R語言中完成線性回歸的方法眾多，本文為了擴(kuò)展至廣義線性回歸分析，介紹glm()完成線性回歸分析[7]，分析代碼如下：

圖1 年齡與脊柱骨密度散點(diǎn)圖

spb<-read.csv("spina_bone.csv") #讀入數(shù)據(jù)

lin_fit <-glm(y～x, family=gaussian(link = "identity"), data=spb) #設(shè)置線性回歸

summary(lin_fit) #輸出結(jié)果

spb$e<-lin_fit $residuals #保存殘差，用于繪制殘差圖

主要分析結(jié)果見表2。

表2 脊柱骨密度研究數(shù)據(jù)線性回歸分析結(jié)果

從分析結(jié)果可知，脊柱骨密度與年齡的數(shù)量關(guān)系為f(xi)=0.494329 + 0.027256xi。但是從殘差圖(圖2)發(fā)現(xiàn)，殘差與年齡存在一定非線性趨勢，表明線性回歸不能完全描述脊柱骨密度與年齡間數(shù)量關(guān)系，導(dǎo)致脊柱骨密度與年齡之間數(shù)量關(guān)系的欠擬合。

圖2 年齡與脊柱骨密度線性回歸殘差圖

3 非參數(shù)回歸及R語言實(shí)現(xiàn)

3.1 非參數(shù)回歸簡介 上述描述發(fā)現(xiàn)脊柱骨密度與年齡之間數(shù)量關(guān)系不是線性，而是復(fù)雜的非線性關(guān)系，即不滿足線性回歸的“線性”要求，需要更為合適的分析方法。將兩者之間的潛在趨勢描述為：

yi=f(xi)+εiE(yi|xi)=f(xi)E(εi)=0

其中f(xi)沒有固定形式，是需要通過數(shù)據(jù)估計(jì)的事先未指定的“光滑”函數(shù)，該分析方法被稱為非參數(shù)回歸分析。常見的非參數(shù)回歸分析方法包含了局部平均、核光滑、局部多項(xiàng)式回歸和樣條光滑等[8-9]。非參數(shù)回歸中結(jié)局變量不局限于高斯分布，可以像參數(shù)回歸一樣擴(kuò)展到指數(shù)分布族。本文介紹基于樣條光滑的非參數(shù)回歸。

進(jìn)行樣條光滑的非參數(shù)回歸需要對解釋變量的變化范圍按節(jié)點(diǎn)分為多個(gè)段，然后進(jìn)行“基變換”[2,8,10]。當(dāng)節(jié)點(diǎn)數(shù)過多時(shí)，導(dǎo)致過擬合現(xiàn)象。過擬合將數(shù)據(jù)的隨機(jī)變異誤以為是系統(tǒng)效應(yīng)，掩蓋真實(shí)的變化趨勢[9]。為了控制非參數(shù)回歸的過擬合現(xiàn)象，需要對光滑函數(shù)進(jìn)行“懲罰”，即在殘差平和中加入“懲罰項(xiàng)”，如下[2,8]：

式中λ(大于等于0)為光滑參數(shù)，通常通過交叉驗(yàn)證(cross-validation，CV)確定。當(dāng)樣本量較大時(shí)，交叉驗(yàn)證耗時(shí)長，可用廣義交叉驗(yàn)證(generalized cross-validation，GCV)替代。

3.2 非參數(shù)回歸的自由度 自由度是非參數(shù)回歸模型復(fù)雜度的評價(jià)指標(biāo)，等同于參數(shù)回歸中自變量的數(shù)量(包含截距)或者自由度-1次多項(xiàng)式回歸，自由度越大，表明模型越復(fù)雜，即解釋變量與結(jié)局變量間關(guān)系越復(fù)雜。如果沒有懲罰項(xiàng)，增加解釋變量“基變換”的節(jié)點(diǎn)，相當(dāng)于增加回歸模型的自由度，會導(dǎo)致過擬合現(xiàn)象。加入“懲罰項(xiàng)”，是保留這些節(jié)點(diǎn)，但是限制影響力，相當(dāng)于減少了回歸模型的自由度。光滑參數(shù)和自由度之間存在單調(diào)的數(shù)學(xué)變換關(guān)系[2]。

3.3 非參數(shù)回歸R語言實(shí)現(xiàn) 非參數(shù)回歸的參數(shù)估計(jì)過程復(fù)雜、計(jì)算量大。在R中進(jìn)行非參數(shù)回歸的函數(shù)主要是gam()函數(shù)[3]，分析代碼如下：

gam_fit<- gam(y～s(x, bs = "cr", k = 22), family=gaussian(link = "identity"), data = spb) #設(shè)置非參數(shù)回歸，內(nèi)部節(jié)點(diǎn)數(shù)量20

summary(gam_fit) #輸出結(jié)果

spb$e<-gam_fit$residuals #保存殘差，用于繪制殘差圖

anova(lin_fit, gam_fit, test = "F") #對非線性進(jìn)行檢驗(yàn)

主要分析結(jié)果見表3。

表3 脊柱骨密度研究數(shù)據(jù)非參數(shù)回歸分析結(jié)果

從分析結(jié)果發(fā)現(xiàn)，截距估計(jì)值為0.929364；光滑函數(shù)對結(jié)局變量有影響，其自由度為5.604，表明年齡與脊柱骨密度間數(shù)量不是線性關(guān)系而是關(guān)系復(fù)雜的非線性關(guān)系；進(jìn)一步對線性參數(shù)回歸和非參數(shù)回歸進(jìn)行比較檢驗(yàn)，發(fā)現(xiàn)兩個(gè)回歸模型存在差異，提示非參數(shù)回歸更優(yōu)(F=23.281，P<0.001)，進(jìn)一步表明年齡與脊柱骨密度間之間不是簡單線性關(guān)系，而是更為復(fù)雜的非線性關(guān)系。從非參數(shù)回歸擬合圖(圖3)，發(fā)現(xiàn)17.5歲之前隨著年齡增長脊柱骨密度增加，而17.5歲之后脊柱骨密度趨于穩(wěn)定；從殘差圖(圖4)，發(fā)現(xiàn)殘差沒有呈現(xiàn)額外趨勢，表明非參數(shù)回歸很好地反映了年齡與脊柱骨密度之間的數(shù)量關(guān)系。

圖3 年齡與脊柱骨密度非參數(shù)回歸擬合結(jié)果

圖4 年齡與脊柱骨密度非參數(shù)回歸殘差圖

4 討論

通過上文闡述，發(fā)現(xiàn)參數(shù)回歸與非參數(shù)回歸有各自的優(yōu)缺點(diǎn)。參數(shù)回歸中函數(shù)部分有固定的表達(dá)式，函數(shù)中參數(shù)的估計(jì)簡單且參數(shù)含義明確，但是參數(shù)回歸函數(shù)形式呆板且對數(shù)據(jù)有著嚴(yán)格要求，難以精確擬合復(fù)雜的非線性關(guān)系，限制了其實(shí)際應(yīng)用。非參數(shù)回歸函數(shù)形式靈活多變沒有固定的表達(dá)式，可以精確擬合復(fù)雜的非線性關(guān)系，但是回歸參數(shù)估計(jì)復(fù)雜且意義模糊。用回歸方法分析具體數(shù)據(jù)時(shí)，不能呆板套用方法，需要對數(shù)據(jù)進(jìn)行精確地“診斷”，才能準(zhǔn)確地分析數(shù)據(jù)，獲得科學(xué)的結(jié)論。筆者推薦如下流程：先繪制散點(diǎn)圖，擬合參數(shù)回歸，擬合非參數(shù)回歸，回歸模型間比較，最后選擇合適回歸分析方法。

非參數(shù)回歸需要控制過擬合現(xiàn)象。先設(shè)置節(jié)點(diǎn)，然后通過光滑參數(shù)或者自由度控制過擬合現(xiàn)象。在分析時(shí)，節(jié)點(diǎn)的數(shù)量和節(jié)點(diǎn)的位置都是需要解決的問題。相對來說，節(jié)點(diǎn)的位置對分析結(jié)果影響較小，為數(shù)據(jù)的等分點(diǎn)即可，節(jié)點(diǎn)的數(shù)量對分析結(jié)果影響較大。多數(shù)情況下，4個(gè)或5個(gè)節(jié)點(diǎn)就能滿足應(yīng)用需求了。大樣本時(shí)，節(jié)點(diǎn)的數(shù)量為20～40[2,9]。相對于節(jié)點(diǎn)的數(shù)量和節(jié)點(diǎn)的位置，光滑參數(shù)與自由度對分析結(jié)果有著決定性的影響，建議用交叉驗(yàn)證或者廣義交叉驗(yàn)證確定。為了滿足更多科研需求，對參數(shù)回歸與非參數(shù)回歸取長補(bǔ)短，在可加模型假設(shè)下，將兩者進(jìn)行結(jié)合，一部分變量是參數(shù)回歸形式，而另一部分是非參數(shù)回歸形式進(jìn)入模型，形成半?yún)?shù)回歸[2-3,8]。半?yún)?shù)回歸是當(dāng)前流行的回歸分析方法。