胡良平

（1.軍事科學(xué)院研究生院，北京 100850；2.世界中醫(yī)藥學(xué)會(huì)聯(lián)合會(huì)臨床科研統(tǒng)計(jì)學(xué)專業(yè)委員會(huì)，北京 100029 *通信作者：胡良平，E-mail：lphu812＠sina.com）

1 概 述［1-2］

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

1.1.1 分位數(shù)

分位數(shù)是一種位置指標(biāo)，一個(gè)特定的分位數(shù)將任何一個(gè)頻數(shù)曲線下的面積（其數(shù)值為1）分為兩部分，若小于等于此分位數(shù)的觀測(cè)值個(gè)數(shù)占全部觀測(cè)值個(gè)數(shù)的比例為1／4，則稱該分位數(shù)為第1四分位數(shù)，記作Q1，同理，還有第2、第3四分位數(shù)，分別記作Q2、Q3；若小于等于此分位數(shù)的觀測(cè)值個(gè)數(shù)占全部觀測(cè)值個(gè)數(shù)的比例為1／10，則稱該分位數(shù)為第1十分位數(shù)，記作D1，同理，還有第2、第3、……、第9十分位數(shù)，分別記作 D2、D3、……、D9；若小于等于此分位數(shù)的觀測(cè)值個(gè)數(shù)占全部觀測(cè)值個(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è)定量因變量的簡(jiǎn)單直線回歸分析中，經(jīng)典統(tǒng)計(jì)學(xué)中需要事先給出如下幾個(gè)假定。

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

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

其三，均值假定。即當(dāng)xi取任何值時(shí)，一律采用y的算術(shù)平均值y-i取代任何一個(gè)yj值，記作y＾=y-i|xi。于是，所有的數(shù)據(jù)點(diǎn)（x，y）可表示為（xi，y＾i）。

其四，獨(dú)立誤差假定。即當(dāng)xi取任何值時(shí)，觀測(cè)值（yj|xi）與其算術(shù)均值 y-i|xi（或 y＾i）之間的偏差為“誤差εi”，假定各誤差εi之間互相獨(dú)立。

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

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

在式（3）和式（4）中，等號(hào)左邊的“y＾i”是在單個(gè)自變量xi（在只有一個(gè)自變量的場(chǎng)合下）或多個(gè)自變量構(gòu)成的向量（Xi）（在有多個(gè)自變量的場(chǎng)合下）取一個(gè)或一組特定數(shù)值的前提條件下，重復(fù)k次觀測(cè)或試驗(yàn)得到了因變量y的k個(gè)觀測(cè)值的“算術(shù)平均值”，簡(jiǎn)稱為“均值”。嚴(yán)格地說，它應(yīng)該被分別表示成如下的形式，見式（5）和式（6）。

對(duì)式（3）而言，y＾i的真實(shí)含義應(yīng)該是式（5）：

對(duì)式（4）而言，y＾i的真實(shí)含義應(yīng)該是式（6）：

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

也就是說，常用的簡(jiǎn)單直線回歸模型或方程和多重線性回歸模型或方程都屬于“均值回歸模型或方程”。“最小二乘原理”可用語(yǔ)言表述如下：

在各xi或Xi條件下，求因變量y的任何一個(gè)觀測(cè)值與其模型對(duì)應(yīng)的估計(jì)值y＾i之間偏差的平方，再將所有各偏差平方求和，稱此“和”為“目標(biāo)函數(shù)Q”，設(shè)法求其最小值，與此同時(shí)，獲得求解模型中參數(shù)估計(jì)值的計(jì)算公式。

值得一提的是：對(duì)基于最小二乘原理構(gòu)造出來的目標(biāo)函數(shù)Q求其最小值時(shí)，需要將Q中的待估參數(shù)視為“變量”，而將各觀測(cè)點(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è)自變量的場(chǎng)合下，中位數(shù)回歸模型見式（7）、中位數(shù)回歸方程見式（8）：

在式（7）和式（8）中，下角標(biāo)中的i代表觀測(cè)的編號(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ù)”所對(duì)應(yīng)的累計(jì)概率，也就是說，“y（0.5）”就是與其對(duì)應(yīng)的“分位數(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ù)對(duì)應(yīng)的回歸模型或回歸方程，見式（9）和式（10）：

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

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

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

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

1.3.1 概述

值得注意的是：在求解“均值回歸模型”中參數(shù)時(shí)，是基于“最小二乘原理”推導(dǎo)出來的正規(guī)方程組并求解得到的；而當(dāng)采用“分位數(shù)”取代“均值”時(shí)，就不適合采取最小二乘原理來構(gòu)建“目標(biāo)函數(shù)Q”，而需要求因變量y的任何一個(gè)觀測(cè)值與其模型對(duì)應(yīng)的估計(jì)值y＾i之間偏差的“絕對(duì)值”，再將所有各點(diǎn)上的偏差絕對(duì)值求和，稱此“和”為“目標(biāo)函數(shù)G”。令式（10）等號(hào)左邊為“W”，則為求解式（10）中參數(shù)估計(jì)值而構(gòu)建的目標(biāo)函數(shù)G如式（11）所示：

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

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

2.1.1 秸稈產(chǎn)量計(jì)算方法 查閱近20篇期刊和論文發(fā)現(xiàn)都采用草谷比計(jì)算法來估算秸稈產(chǎn)量。其計(jì)算方法為：

由 SAS軟件的幫助信息［從 SAS／STAT的QUANTREG過程的“details”（即詳細(xì)情況）菜單進(jìn)入，其第 2行為“Optimization Algorithms”（優(yōu)化算法）］可知，求式（11）最小值的計(jì)算方法有三種，分別為“簡(jiǎn)單算法”“內(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）的觀測(cè)值。見表 1。

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

續(xù)表1：

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

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

比較式（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ù)中，對(duì)因變量y而言，后50個(gè)y值明顯大于前950個(gè)y值，故屬于“異常值”，它們所對(duì)應(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 對(duì)資料中因變量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 對(duì)表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 對(duì)資料進(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語(yǔ)句中的選項(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 Function

Quantile 0.25

Objective Function 1282.8452

Predicted Value at Mean 9.6857

Parameter Estimates

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

Quantile and Objective Function

Parameter Estimates

以上是“τ=0.50”時(shí)對(duì)應(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 Function

Parameter Estimates

以上呈現(xiàn)的是“τ=0.75”時(shí)對(duì)應(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|Intercept 1 14.48953 0.62584 23.15 ＜.0001 x1 1 4.39030 0.62997 6.97 ＜.0001 x2 1 2.50293 0.60204 4.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ì)如下：

LTSParameter Estimates

以上結(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）：

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