胡良平

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

1 概 述[1-2]

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

1.1.1 分位數(shù)

分位數(shù)是一種位置指標(biāo)，一個特定的分位數(shù)將任何一個頻數(shù)曲線下的面積(其數(shù)值為1)分為兩部分，若小于等于此分位數(shù)的觀測值個數(shù)占全部觀測值個數(shù)的比例為1/4，則稱該分位數(shù)為第1四分位數(shù)，記作Q1，同理，還有第2、第3四分位數(shù)，分別記作Q2、Q3；若小于等于此分位數(shù)的觀測值個數(shù)占全部觀測值個數(shù)的比例為1/10，則稱該分位數(shù)為第1十分位數(shù)，記作D1，同理，還有第2、第3、……、第9十分位數(shù)，分別記作D2、D3、……、D9； 若小于等于此分位數(shù)的觀測值個數(shù)占全部觀測值個數(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ù)。通過給出一組資料的若干個分位數(shù)，可初步描述該組資料的離散程度和分布概況，故在實際工作中，常用百分位數(shù)法確定服從偏態(tài)分布資料的醫(yī)學(xué)指標(biāo)的正常值范圍。

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

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

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

其二，同方差假定。即當(dāng)xi取任何值時，其對應(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ù)并按最小二乘原理估計出其參數(shù)的數(shù)值，然后忽略模型中的誤差項，就可得到“均值回歸方程”，見式(3)和式(4)：

(3)

(4)

(5)

(6)

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

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

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

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

在前面構(gòu)建“均值回歸模型或方程”的“均值假定”中，若用“中位數(shù)”取代“算術(shù)平均值”，就得到了“中位數(shù)回歸模型或回歸方程”，在有m個自變量的場合下，中位數(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代表觀測的編號，i=1,2,……,n；“0”代表“常數(shù)項”；“1、2、……、m”分別代表“第1個”“第2個”……“第m個”自變量；“0.5”代表“第50百分位數(shù)”或“第2四分位數(shù)”或“中位數(shù)”所對應(yīng)的累計概率，也就是說，“y(0.5)”就是與其對應(yīng)的“分位數(shù)”，即“中位數(shù)”。

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

由前面關(guān)于“分位數(shù)”的概念可知：“中位數(shù)”是一個特定的“分位數(shù)”。若將上面的“中位數(shù)回歸模型或方程”中的“中位數(shù)”替換成任意一個分位數(shù)“yτ|Xi”，則所得到的回歸模型就被稱為“分位數(shù)回歸模型或方程”了。前述的“yτ|Xi”的含義是：在自變量x取特定值xi的條件下，求出全部因變量y的第τ分位數(shù)，0<τ<1。當(dāng)τ=0.5時，就是“第0.5分位數(shù)”或“中位數(shù)”；當(dāng)τ=0.25時，就是“第0.25分位數(shù)”或“第1/4分位數(shù)”，也叫做“第1四分位數(shù)”；當(dāng)τ=0.75時，就是“第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ù)回歸分析中，每給定一個“τ值”，就可求出一個相應(yīng)的“分位數(shù)回歸模型或方程”。故對于一個給定的資料來說，因τ在開區(qū)間(0，1)范圍內(nèi)有無窮多個取值，就有無窮多個“分位數(shù)回歸模型或方程”。其中，最常見的分位數(shù)回歸模型或方程為“第1四分位數(shù)回歸模型或方程”“第2四分位數(shù)回歸模型(即中位數(shù)回歸模型)或方程”和“第3四分位數(shù)回歸模型或方程”。

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

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

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

1.3.1 概述

(11)

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

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

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

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

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

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

【例1】假定有一個總樣本含量n=1 000的數(shù)據(jù)集中包含5%異常點的資料，每組數(shù)據(jù)(即每個個體)包含三個變量(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省略編號為11～990的980行數(shù)據(jù)；在全部1000行數(shù)據(jù)中，最后50行數(shù)據(jù)為“異常點”，占5%

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

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

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

(12)

y=100+e

(13)

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

【問題】試擬合表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)性檢驗并繪制其頻數(shù)直方圖，直觀了解y的頻數(shù)分布情況。

proc univariate data=a;

var y;

histogram y;

run;

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

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

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

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

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

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

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

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

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

run;

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

【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ù)估計及檢驗結(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”時對應(yīng)的“第2/4四分位數(shù)回歸模型”(即中位數(shù)回歸模型)的參數(shù)估計及檢驗結(jié)果，結(jié)果表明：截距項=10.0364、兩個自變量的斜率分別為5.0106和3.0294，參考它們的理論值[見式(12)]分別為10、5和3，說明在因變量不服從正態(tài)分布且資料中含有5%異常點時，采用“中位數(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”時對應(yīng)的“第3/4四分位數(shù)回歸模型”的參數(shù)估計及檢驗結(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)計學(xué)意義；但復(fù)相關(guān)系數(shù)的平方R2=0.0634非常小，表明此二重線性模型的實用價值并不高；模型中三個參數(shù)的估計與假設(shè)檢驗結(jié)果如下：

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

以上結(jié)果表明：截距項=14.48953，兩個自變量的斜率分別為4.39030、2.50293，參考它們的理論值[見式(12)]分別為10、5和3，說明在因變量不服從正態(tài)分布且資料中含有5%異常點時，直接基于普通最小二乘原理創(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ù)估計如下：

LTS Parameter EstimatesParameterDFEstimateIntercept110.0054x115.0240x213.0598

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

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

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

(14)(OLS估計法)

(15)(LTS估計法)

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

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