分位數(shù)模型回歸分析

胡良平

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

1 概 述[1-2]

1.1 分位數(shù)回歸模型

1.1.1 分位數(shù)

分位數(shù)是一種位置指標(biāo)，一個(gè)特定的分位數(shù)將任何一個(gè)頻數(shù)曲線下的面積(其數(shù)值為1)分為兩部分，若小于等于此分位數(shù)的觀測值個(gè)數(shù)占全部觀測值個(gè)數(shù)的比例為1/4，則稱該分位數(shù)為第1四分位數(shù)，記作Q1，同理，還有第2、第3四分位數(shù)，分別記作Q2、Q3；若小于等于此分位數(shù)的觀測值個(gè)數(shù)占全部觀測值個(gè)數(shù)的比例為1/10，則稱該分位數(shù)為第1十分位數(shù)，記作D1，同理，還有第2、第3、……、第9十分位數(shù)，分別記作D2、D3、……、D9； 若小于等于此分位數(shù)的觀測值個(gè)數(shù)占全部觀測值個(gè)數(shù)的比例為1/100，則稱該分位數(shù)為第1百分位數(shù)，記作P1，同理，還有第2、第3、……、第99百分位數(shù)，分別記作P2、P3、……、P99。

顯然，第1四分位數(shù)=第25百分位數(shù)，即Q1=P25；第2四分位數(shù)=第5十分位數(shù)=第50百分位數(shù)=中位數(shù)，即Q2=D5=P50=M(代表中位數(shù))；第3四分位數(shù)=第75百分位數(shù)，即Q3=P75。如此，常用百分位數(shù)代替四分位數(shù)和十分位數(shù)。通過給出一組資料的若干個(gè)分位數(shù)，可初步描述該組資料的離散程度和分布概況，故在實(shí)際工作中，常用百分位數(shù)法確定服從偏態(tài)分布資料的醫(yī)學(xué)指標(biāo)的正常值范圍。

1.1.2 均值回歸模型與均值回歸方程

在僅有一個(gè)定量自變量和一個(gè)定量因變量的簡單直線回歸分析中，經(jīng)典統(tǒng)計(jì)學(xué)中需要事先給出如下幾個(gè)假定。

其一，正態(tài)假定。即當(dāng)自變量x在其取值區(qū)間內(nèi)取定任何一個(gè)特定值xi時(shí)，因變量y可以有一組取值與其對應(yīng)。例如很多身高x為165cm的成年人其體重y是不等的。如此多的y值就會(huì)有一個(gè)概率分布，粗略地說，該分布可能是正偏態(tài)的、對稱的(最理想的為正態(tài)的)或負(fù)偏態(tài)的，簡單記作P(y|xi)，經(jīng)典統(tǒng)計(jì)學(xué)假定P(y|xi)為“正態(tài)分布”，對所有的xi都成立。

其二，同方差假定。即當(dāng)xi取任何值時(shí)，其對應(yīng)的上述正態(tài)分布P(y|xi)的方差均相等，稱為“同方差”或“方差齊性”。否則，就稱為“異方差”或“方差不齊”。

基于以上假定且按“普通最小二乘原理或最小平方原理”[3]構(gòu)建的簡單直線回歸模型和多重線性回歸模型都被稱為“均值回歸模型”，分別見式(1)和式(2)：

yi=β0+β1xi+εii=1,2,…,n

(1)

yi=β0+β1x1i+β2x2i+…+βmxmi+εii=1,2,…,n

(2)

在模型(1)和(2)中，基于樣本數(shù)據(jù)并按最小二乘原理估計(jì)出其參數(shù)的數(shù)值，然后忽略模型中的誤差項(xiàng)，就可得到“均值回歸方程”，見式(3)和式(4)：

(3)

(4)

(5)

(6)

在式(6)中，大寫的字母“X”代表由多個(gè)自變量構(gòu)成的向量。

也就是說，常用的簡單直線回歸模型或方程和多重線性回歸模型或方程都屬于“均值回歸模型或方程”?！白钚《嗽怼笨捎谜Z言表述如下：

值得一提的是：對基于最小二乘原理構(gòu)造出來的目標(biāo)函數(shù)Q求其最小值時(shí)，需要將Q中的待估參數(shù)視為“變量”，而將各觀測點(diǎn)上變量的取值(xi,yi)(一個(gè)自變量的情況下)或(x1i,x2i,…,xmi,yi)(m個(gè)自變量的情況下)視為“常數(shù)”，利用高等數(shù)學(xué)中求函數(shù)“極大值”和“極小值”的方法，即求函數(shù)Q關(guān)于未知參數(shù)的偏導(dǎo)數(shù)等多個(gè)計(jì)算步驟來實(shí)現(xiàn)目的。事實(shí)上，求得的偏導(dǎo)數(shù)又是關(guān)于研究問題中“變量(x，y)或(x1,x2,…,xm,y)”的函數(shù)，故稱其為“導(dǎo)函數(shù)”。令求得的導(dǎo)函數(shù)為零，就得到若干個(gè)方程，其個(gè)數(shù)為所構(gòu)建的回歸模型中參數(shù)的個(gè)數(shù)。在直線回歸模型中，有兩個(gè)參數(shù)，即截距和斜率；在含有m個(gè)自變量的多重線性回歸模型中，有(m+1)個(gè)參數(shù)。將這些線性方程聯(lián)立在一起，就構(gòu)成了“線性方程組”。求出該方程組的解，就得到了求解待估計(jì)參數(shù)的計(jì)算公式。

1.1.3 中位數(shù)回歸模型與中位數(shù)回歸方程

在前面構(gòu)建“均值回歸模型或方程”的“均值假定”中，若用“中位數(shù)”取代“算術(shù)平均值”，就得到了“中位數(shù)回歸模型或回歸方程”，在有m個(gè)自變量的場合下，中位數(shù)回歸模型見式(7)、中位數(shù)回歸方程見式(8)：

y(0.5)i|(X=Xi)=β0(0.5)+β1(0.5)ix1i+β2(0.5)ix2i+…+βm(0.5)ixmi+ε(0.5)i

(7)

(8)

在式(7)和式(8)中，下角標(biāo)中的i代表觀測的編號(hào)，i=1,2,……,n；“0”代表“常數(shù)項(xiàng)”；“1、2、……、m”分別代表“第1個(gè)”“第2個(gè)”……“第m個(gè)”自變量；“0.5”代表“第50百分位數(shù)”或“第2四分位數(shù)”或“中位數(shù)”所對應(yīng)的累計(jì)概率，也就是說，“y(0.5)”就是與其對應(yīng)的“分位數(shù)”，即“中位數(shù)”。

1.1.4 分位數(shù)回歸模型與分位數(shù)回歸方程

由前面關(guān)于“分位數(shù)”的概念可知：“中位數(shù)”是一個(gè)特定的“分位數(shù)”。若將上面的“中位數(shù)回歸模型或方程”中的“中位數(shù)”替換成任意一個(gè)分位數(shù)“yτ|Xi”，則所得到的回歸模型就被稱為“分位數(shù)回歸模型或方程”了。前述的“yτ|Xi”的含義是：在自變量x取特定值xi的條件下，求出全部因變量y的第τ分位數(shù)，0<τ<1。當(dāng)τ=0.5時(shí)，就是“第0.5分位數(shù)”或“中位數(shù)”；當(dāng)τ=0.25時(shí)，就是“第0.25分位數(shù)”或“第1/4分位數(shù)”，也叫做“第1四分位數(shù)”；當(dāng)τ=0.75時(shí)，就是“第0.75分位數(shù)”或“第3/4分位數(shù)”，也叫做“第3四分位數(shù)”。

在式(7)和式(8)中，若將“0.5”改換成“τ∈(0，1)”(其含義是0<τ<1)，就得到了與第τ分位數(shù)對應(yīng)的回歸模型或回歸方程，見式(9)和式(10)：

y(τ)i|(X=Xi)=β0(τ)+β1(τ)ix1i+β2(τ)ix2i+…+βm(τ)ixmi+ε(τ)i

(9)

(10)

在分位數(shù)回歸分析中，每給定一個(gè)“τ值”，就可求出一個(gè)相應(yīng)的“分位數(shù)回歸模型或方程”。故對于一個(gè)給定的資料來說，因τ在開區(qū)間(0，1)范圍內(nèi)有無窮多個(gè)取值，就有無窮多個(gè)“分位數(shù)回歸模型或方程”。其中，最常見的分位數(shù)回歸模型或方程為“第1四分位數(shù)回歸模型或方程”“第2四分位數(shù)回歸模型(即中位數(shù)回歸模型)或方程”和“第3四分位數(shù)回歸模型或方程”。

1.2 分位數(shù)模型回歸分析應(yīng)用的場合

當(dāng)擬做多重線性回歸分析的原始數(shù)據(jù)中的定量因變量不服從正態(tài)分布、有時(shí)還存在異方差、資料中存在一定比例的“異常點(diǎn)”且自變量間不存在嚴(yán)重多重共線性時(shí)，采用此方法構(gòu)建多重線性回歸模型，可以最大限度地消除資料中違反經(jīng)典回歸分析的“部分或全部假定”對建模結(jié)果造成的影響。

1.3 分位數(shù)模型回歸分析的計(jì)算原理

1.3.1 概述

(11)

設(shè)法求出式(11)的最小值，與此同時(shí)，獲得模型中參數(shù)的估計(jì)值，這樣求出的參數(shù)估計(jì)值被稱為“分位數(shù)回歸參數(shù)估計(jì)值”。

1.3.2 式(11)最小值的計(jì)算方法

由SAS軟件的幫助信息[從SAS/STAT的QUANTREG過程的“details”(即詳細(xì)情況)菜單進(jìn)入，其第2行為“Optimization Algorithms”(優(yōu)化算法)]可知，求式(11)最小值的計(jì)算方法有三種，分別為“簡單算法”“內(nèi)部點(diǎn)算法”和“光滑算法”，在樣本含量n<5 000且自變量個(gè)數(shù)m<50時(shí)，三種算法的參數(shù)估計(jì)結(jié)果是基本相同的。

事實(shí)上，上述提及的三種算法都相當(dāng)復(fù)雜，涉及到復(fù)雜的矩陣運(yùn)算和反復(fù)迭代計(jì)算過程。換言之，由式(11)無法給出與各種分位數(shù)回歸模型或方程中各參數(shù)的解析式，即計(jì)算公式。故具體分析時(shí)，建議采用統(tǒng)計(jì)軟件來實(shí)現(xiàn)計(jì)算。因篇幅所限，詳細(xì)的計(jì)算方法此處不再贅述。

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

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

【例1】假定有一個(gè)總樣本含量n=1 000的數(shù)據(jù)集中包含5%異常點(diǎn)的資料，每組數(shù)據(jù)(即每個(gè)個(gè)體)包含三個(gè)變量(x1,x2,y)的觀測值。見表1。

表1 某資料中首尾各10組數(shù)據(jù)(n=1000)

續(xù)表1：

8-0.481391.3382111.2089980.55630-1.52099105.3409-0.23384-0.859546.0609990.051370.2431995.097100.00900-0.243769.6611000-1.539940.96521105.054

注：表1省略編號(hào)為11～990的980行數(shù)據(jù)；在全部1000行數(shù)據(jù)中，最后50行數(shù)據(jù)為“異常點(diǎn)”，占5%

【特別說明】例1是人為構(gòu)造的，它來自SAS 9.3的QUANTREG過程中的“樣例”。三個(gè)變量“x1、x2和y”沒有實(shí)際的專業(yè)含義，僅為了造出一個(gè)樣本含量為1 000且含5%異常點(diǎn)的數(shù)據(jù)集。

表1中數(shù)據(jù)構(gòu)造的方法如下：設(shè)定x1和x2及測量誤差e都是服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)變量(其均值為0、方差為1)，前950個(gè)y的數(shù)值按下面的模型(12)計(jì)算出來；后50個(gè)y的數(shù)值按下面的模型(13)計(jì)算出來：

y=10+5*x1+3*x2+0.5 *e

(12)

y=100+e

(13)

比較式(12)與式(13)可知：y的前950個(gè)數(shù)據(jù)中的每一個(gè)都在基數(shù)“10”的基礎(chǔ)上再加上三項(xiàng)并不大的數(shù)值，其平均值約為“10+5+3=18”；而y的后50個(gè)數(shù)據(jù)中的每一個(gè)都在基數(shù)“100”的基礎(chǔ)上再加上一個(gè)隨機(jī)誤差，其平均值約為100。由此可知：表1的1000行數(shù)據(jù)中，對因變量y而言，后50個(gè)y值明顯大于前950個(gè)y值，故屬于“異常值”，它們所對應(yīng)的那50行數(shù)據(jù)點(diǎn)就屬于“異常點(diǎn)”了。

【問題】試擬合表1中y依賴x1、x2變化的二重線性回歸模型。

2.2 所需要的SAS程序

2.2.1 產(chǎn)生數(shù)據(jù)集

先用下面的一段SAS數(shù)據(jù)步程序產(chǎn)生表1中的1000行3列數(shù)據(jù)，創(chuàng)建數(shù)據(jù)集a。

data a (drop=i);

do i=1 to 1000;

x1=rannor(1234);

x2=rannor(1234);

e=rannor(1234);

if i>950 then y=100+10*e;

else y=10+5*x1+3*x2+0.5*e;

output;

end;

run;

/* 以上程序產(chǎn)生1000行數(shù)據(jù)(x1,x2,y)，其中，有5%的是異常值 */

以上程序產(chǎn)生的數(shù)據(jù)見表1。

2.2.2 對資料中因變量y的分布情況進(jìn)行探索性分析

采用下面的SAS程序?qū)σ蜃兞縴進(jìn)行正態(tài)性檢驗(yàn)并繪制其頻數(shù)直方圖，直觀了解y的頻數(shù)分布情況。

proc univariate data=a;

var y;

histogram y;

run;

以上SAS程序輸出結(jié)果見表2、圖1。

表2 對表1中1000個(gè)y值是否服從正態(tài)分布檢驗(yàn)結(jié)果

由表2第1行“W檢驗(yàn)”可知，本例中y值不服從正態(tài)分布。

圖1顯示：1 000個(gè)y值中的950個(gè)在圖中左邊部分，呈正偏態(tài)分布；中間出現(xiàn)了較大一段空缺，還有50個(gè)y值位于圖中右邊部分，其數(shù)值波動(dòng)約為72～120。

2.2.3 對資料進(jìn)行分位數(shù)模型回歸分析

下面基于“τ=0.25、0.50和0.75”分別用“分位數(shù)回歸分析法”來創(chuàng)建二重線性回歸模型：

圖1 表1中全部1000個(gè)y值的頻數(shù)分布直方圖

proc quantreg algorithm=smooth(ratio=.5) data=a;

model y = x1 x2/quantile=0.25 0.50 0.75;

run;

【SAS程序說明】model語句中的選項(xiàng)“quantile=0.25 0.50 0.75”是要求SAS系統(tǒng)分別創(chuàng)建分位數(shù)為0.25、0.50(中位數(shù))和0.75的3個(gè)二重線性回歸模型。

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

Quantile and Objective FunctionQuantile0.25Objective Function1282.8452Predicted Value at Mean9.6857

Parameter EstimatesParameterDFEstimateStandard Error95% Confidence Limitst ValuePr > |t|Intercept19.69570.01989.65699.7345490.27<.0001x115.00830.01954.97015.0466257.09<.0001x213.01700.01482.98803.0460204.22<.0001

以上是“τ=0.25”對應(yīng)的“第1/4四分位數(shù)回歸模型”的參數(shù)估計(jì)及檢驗(yàn)結(jié)果

Quantile and Objective FunctionQuantile0.5Objective Function2441.1927Predicted Value at Mean10.0259

Parameter EstimatesParameterDFEstimateStandard Error95% Confidence Limitst ValuePr > |t|Intercept110.03640.02199.993410.0794457.96<.0001x115.01060.01944.97255.0487257.90<.0001x213.02940.01862.99303.0658163.31<.0001

以上是“τ=0.50”時(shí)對應(yīng)的“第2/4四分位數(shù)回歸模型”(即中位數(shù)回歸模型)的參數(shù)估計(jì)及檢驗(yàn)結(jié)果，結(jié)果表明：截距項(xiàng)=10.0364、兩個(gè)自變量的斜率分別為5.0106和3.0294，參考它們的理論值[見式(12)]分別為10、5和3，說明在因變量不服從正態(tài)分布且資料中含有5%異常點(diǎn)時(shí)，采用“中位數(shù)回歸分析”創(chuàng)建的二重線性回歸模型與其“原模型”之間的吻合程度非常高。

Quantile and Objective FunctionQuantile0.75Objective Function3512.2464Predicted Value at Mean10.4106

Parameter EstimatesParameterDFEstimateStandard Error95% Confidence Limitst ValuePr > |t|Intercept110.42050.024710.372010.4690421.73<.0001x114.99930.02754.94535.0532181.90<.0001x213.00930.02802.95443.0641107.66<.0001

以上呈現(xiàn)的是“τ=0.75”時(shí)對應(yīng)的“第3/4四分位數(shù)回歸模型”的參數(shù)估計(jì)及檢驗(yàn)結(jié)果。

2.3 小結(jié)[3-4]

2.3.1 采用普通最小二乘法構(gòu)建二重線性回歸模型

proc reg data=a;

model y=x1 x2;

run;

以上SAS程序的主要輸出結(jié)果如下：

總模型的F=33.76，P<0.0001，表明所創(chuàng)建的二重線性回歸模型具有統(tǒng)計(jì)學(xué)意義；但復(fù)相關(guān)系數(shù)的平方R2=0.0634非常小，表明此二重線性模型的實(shí)用價(jià)值并不高；模型中三個(gè)參數(shù)的估計(jì)與假設(shè)檢驗(yàn)結(jié)果如下：

變量自由度參數(shù)估計(jì)值標(biāo)準(zhǔn)誤差t 值Pr > |t|Intercept114.489530.6258423.15<.0001x114.390300.629976.97<.0001x212.502930.602044.16<.0001

以上結(jié)果表明：截距項(xiàng)=14.48953，兩個(gè)自變量的斜率分別為4.39030、2.50293，參考它們的理論值[見式(12)]分別為10、5和3，說明在因變量不服從正態(tài)分布且資料中含有5%異常點(diǎn)時(shí)，直接基于普通最小二乘原理創(chuàng)建的二重線性回歸模型很不理想。

2.3.2采用“穩(wěn)健回歸分析”創(chuàng)建二重線性回歸模型

proc robustreg data=a method=lts seed=100;

model y=x1 x2;

run;

以上SAS程序的主要輸出結(jié)果如下：

復(fù)相關(guān)系數(shù)的平方R2=0.9933非常大，表明此二重線性回歸模型具有很高的實(shí)用價(jià)值，其參數(shù)估計(jì)如下：

LTS Parameter EstimatesParameterDFEstimateIntercept110.0054x115.0240x213.0598

以上結(jié)果表明：截距項(xiàng)=10.0054，兩個(gè)自變量的斜率分別為5.0240和3.0598，參考它們的理論值[見式(12)]分別為10、5和3，說明在因變量不服從正態(tài)分布且資料中含有5%異常點(diǎn)時(shí)，基于最小截平方和(Least trimmed squares,LTS)法的“穩(wěn)健回歸分析”創(chuàng)建的二重線性回歸模型與其“原模型”非常吻合。

2.3.3 總結(jié)三種建模方法所得結(jié)果的差異

在因變量不服從正態(tài)分布且資料中存在異常點(diǎn)時(shí)，采用SAS中REG過程并基于普通最小二乘原理直接創(chuàng)建多重線性回歸模型的效果很差，而基于“穩(wěn)健回歸分析法”和“分位數(shù)回歸分析法”得到的結(jié)果非常接近。事實(shí)上，本例中的數(shù)據(jù)是基于上面的二重線性回歸模型(12)產(chǎn)生的。這個(gè)模型意味著：截距項(xiàng)為10、x1前的回歸系數(shù)為5、x2前的回歸系數(shù)為3，在此基礎(chǔ)上，加一個(gè)隨機(jī)誤差的二分之一。此處的“隨機(jī)誤差”是服從均值為0、方差為1的正態(tài)分布的隨機(jī)誤差?；谏厦娴挠?jì)算結(jié)果，可以寫出三個(gè)二重線性回歸模型的具體表達(dá)式如下，見式(14)、式(15)、式(16)：

(14)(OLS估計(jì)法)

(15)(LTS估計(jì)法)

(16)(分位數(shù)估計(jì)法)

【結(jié)論】以模型(12)為“金標(biāo)準(zhǔn)”，模型(14)偏離很遠(yuǎn)；模型(15)和(16)的質(zhì)量都很高。相比之下，中位數(shù)回歸分析的結(jié)果稍優(yōu)于穩(wěn)健回歸分析法(以“回歸系數(shù)與其真值之間的偏差最小”為評價(jià)依據(jù))。