何鳳霞, 王鳳竹

?

分位數(shù)回歸及其在R中的實(shí)現(xiàn)

何鳳霞*, 王鳳竹

(華北電力大學(xué) 數(shù)理學(xué)院, 北京, 102206)

針對R中沒有函數(shù)可以實(shí)現(xiàn)模型似然比檢驗(yàn)、模型的擬合度等問題,編寫了擬似然檢驗(yàn)、擬合度的R代碼, 運(yùn)用實(shí)例給出分位數(shù)回歸在R中實(shí)現(xiàn)的一般流程, 為以后分位數(shù)回歸的算法研究提供了重要參考.

分位數(shù)回歸; R; 擬似然比; 擬合度

經(jīng)典回歸與分位數(shù)回歸的主要區(qū)別是: 經(jīng)典回歸描述自變量對因變量的條件均值影響, 而分位數(shù)回歸描述自變量對于整個(gè)因變量的條件變化影響; 經(jīng)典回歸的隨機(jī)項(xiàng)需來自均值為0且同方差的正態(tài)分布, 而分位數(shù)回歸對于隨機(jī)項(xiàng)不需具體分布的假設(shè); 經(jīng)典回歸易受離群點(diǎn)的影響, 而分位數(shù)回歸對于離群點(diǎn)不敏感等.

1978年, Koenker和Bassett[1]在最小絕對偏差估計(jì)理論的基礎(chǔ)上首次提出了分位數(shù)回歸的概念, 分位數(shù)回歸在近20年廣泛應(yīng)用于經(jīng)濟(jì)、醫(yī)學(xué)、生物等領(lǐng)域. 2006年, 李育安[2]介紹了分位數(shù)回歸的概念及其應(yīng)用; 荀鵬程等[3]將其應(yīng)用于預(yù)測SARS的發(fā)病率.

目前，能實(shí)現(xiàn)分位數(shù)回歸的軟件主要有SAS、Eviews、Stata、R. R軟件是一種開源、免費(fèi)的優(yōu)秀統(tǒng)計(jì)軟件, 并且可以和全球一流統(tǒng)計(jì)專家討論; 很多統(tǒng)計(jì)學(xué)的前沿方法均能實(shí)現(xiàn), 其中quantreg包可以實(shí)現(xiàn)分位數(shù)回歸.

本文主要介紹分位數(shù)回歸(以線性為例)的概念、參數(shù)估計(jì)、參數(shù)檢驗(yàn)以及模型的檢驗(yàn); 其次給出分位數(shù)回歸參數(shù)估計(jì)、參數(shù)檢驗(yàn)用R軟件實(shí)現(xiàn)的函數(shù), 對于R中沒有函數(shù)可以實(shí)現(xiàn)模型似然比檢驗(yàn)、模型的擬合度給出自編R代碼; 最后以R自帶數(shù)據(jù)為例, 給出一般分位數(shù)回歸在R中實(shí)現(xiàn)的步驟.

1 分位數(shù)回歸的簡介

2 參數(shù)的估計(jì)

2.1 單純形算法

單純形算法[2]是由Koenker提出, 該算法適合樣本量不大和自變量個(gè)數(shù)不多的變量, 當(dāng)數(shù)據(jù)中存在大量離群點(diǎn)時(shí), 單純形算法估計(jì)出來的參數(shù)穩(wěn)定性比較好, 但是在處理大量數(shù)據(jù)時(shí)運(yùn)算的速度會顯著降低[3]. R中實(shí)現(xiàn)函數(shù)為rq(formula, tau = 0.5, method = "br"), 其中formula表示公式對象; tau為分位點(diǎn), 默認(rèn)為中位數(shù)回歸(即0.5分位數(shù)回歸), 如建立自變量對因變量的0.9分位數(shù)回歸, rq(formula, tau = 0.9); method表示估計(jì)參數(shù)的方法, 單純形算法用"br"表示, 其實(shí)R中默認(rèn)算法是單純形算法, 一般樣本容量不超過5 000、自變量個(gè)數(shù)不大于20, 經(jīng)常用單純形算法.

2.2 內(nèi)點(diǎn)法

單純形算法處理大樣本時(shí)效率比較低, Barrodale、Koenker提出了內(nèi)點(diǎn)法, 適合樣本量較大, 且自變量個(gè)數(shù)少的樣本數(shù)據(jù)[4]. 關(guān)于單純形算法與內(nèi)點(diǎn)法的比較, Koenker模擬數(shù)據(jù)[5]比較了自變量個(gè)數(shù)為= 4, 8, 16, 樣本容量= 200, 400, 800, 1 200, 2 000, 4 000, 8 000, 12 000的情形; 得到當(dāng)樣本量比較少時(shí), 單純形算法優(yōu)于內(nèi)點(diǎn)法, 隨著樣本量的增加內(nèi)點(diǎn)法優(yōu)于單純形算法. 在R中內(nèi)點(diǎn)法的實(shí)現(xiàn)函數(shù)為rq (formula, tau, method = "fn").

除了單純形算法和內(nèi)點(diǎn)法在R中可以實(shí)現(xiàn)外, R也可以實(shí)現(xiàn)Portnoy、Koenker提出的預(yù)處理內(nèi)點(diǎn)法等[4], 命令分別為rq(~, method = "pfn").

3 參數(shù)檢驗(yàn)

本文主要介紹參數(shù)檢驗(yàn)常用方法的3種: 誤差為獨(dú)立同分布、誤差為非獨(dú)立同分布、Bootstrap法.

3.1 誤差為獨(dú)立同分布

R中的命令為summary.rq(object, se = "iid", hs = T), 其中object是由上面命令rq返回的對象; se表示參數(shù)檢驗(yàn)使用的方法; 當(dāng)hs = T(默認(rèn))使用Hall-Sheater寬帶, hs = F時(shí), 使用Boginger寬帶, 需要注意的是R中稀疏函數(shù)的估計(jì)使用RogerKoenker改進(jìn)了差分方法.

3.2 誤差為非獨(dú)立同分布

對于誤差為非獨(dú)立同分布時(shí), 參數(shù)的極限分布比較復(fù)雜, 如位移尺度假設(shè)模型:

3.3 Bootstrap法

Bootstrap法通過重新抽樣來估計(jì)參數(shù), 常用的方法有-Bootstrap、馬爾可夫鏈邊際Bootstrap.

3.3.1-成對Bootstrap

-Bootstrap是Bootstrap最一般的形式, 它首先對(,)進(jìn)行重新抽樣(抽取的樣本數(shù)小于等于初始的樣本數(shù)); 再計(jì)算分位數(shù)回歸的系數(shù)估計(jì)值, 重復(fù)進(jìn)行次抽樣得到個(gè)系數(shù)估計(jì)值; 最后得到參數(shù)的漸近分布協(xié)方差. 在R中通過summary(object, se = "boot", bsmethod = "xy")實(shí)現(xiàn).

3.32 馬爾可夫鏈邊際Bootstrap

He和Hu提出馬爾可夫鏈邊際Bootstrap[6], 它將多維線性規(guī)劃問題轉(zhuǎn)化為一維, 這些一維的解構(gòu)成馬爾可夫鏈, 簡化運(yùn)算. R中實(shí)現(xiàn)函數(shù)為summary(object, se = "boot", bsmethod = "mcmcb"). R中除了能實(shí)現(xiàn)上面2種方法外, 還可以實(shí)現(xiàn)Parzen和Wei and Ying提出的pwy法[7]等.

4 模型檢驗(yàn)及擬合度

4.1 模型的檢驗(yàn)

對于模型的檢驗(yàn), Koenker提出擬似然比(quasi-likelood-ration)[8]檢驗(yàn)整個(gè)模型的顯著性. 它假設(shè)所有參數(shù)均為0. 構(gòu)造的統(tǒng)計(jì)量為:

4.2 模型的擬合度

分位數(shù)回歸的擬合度[8]是Koerker、Machado提出的, 它的值介于0與1之間, 其表達(dá)式為:

5 實(shí)例

以R中quantreg包自帶數(shù)據(jù)engel[9—11](這個(gè)數(shù)據(jù)是Roenker給出的, 它包含收入income、食物支出foodexp變量, 253個(gè)觀察值)為例, 給出R中分位數(shù)回歸的一般步驟.

第1步: 參數(shù)估計(jì).

采用單純形法求解參數(shù)估計(jì)tau = 0.5的回歸系數(shù). 命令:

>library(quantreg)#調(diào)用quantreg包

>data(engel)#讀取engel數(shù)據(jù)

>attach(engel)#連接數(shù)據(jù)框的變量名到內(nèi)存中, 以便直接調(diào)用變量

>fit=rq(foodexp~income, tau = 0.5, method = "br")#用單純形法估計(jì)中位數(shù)回歸的參數(shù)

>fit#查看回歸估計(jì)參數(shù)

Coefficients:

Intercept income

81.482 0.560

結(jié)果解釋: 中位數(shù)回歸的常數(shù)項(xiàng)為81.482, 系數(shù)為0.56.

第2步: 參數(shù)的檢驗(yàn).

采用bootstrap的X-Y方法. R中的命令:

>summary(fit, se = "boot", bsmethod = "xy")

Value Std. Error t value Pr(> |t|)

Intercept 81.482 6.932 3.025 0.003

Income0.560 0.034 16.579 0.000

結(jié)果解釋: 可以得到常數(shù)項(xiàng)的= 0.003 < 0.05, 系數(shù)的= 0.000 < 0.05, 所以在顯著水平為0.05下參數(shù)均顯著;

第3步: 模型的檢驗(yàn).

>source("C:/QLR.r")#調(diào)用QLR.r

>QLR(fit, hs = T, mod = F)

$tau

[1] 0.5

$QLR.statistic

[1] 437.4879

$p.value

[1] 0

結(jié)果解釋: tau = 0.5的擬似然比統(tǒng)計(jì)量值為437.4879,值為0. 由于值小于0.05, 所以在顯著水平0.05下顯著.

第4步: 擬合度的計(jì)算.

>source("C:/R1.r")#調(diào)用R1.r

>R1(fit)#計(jì)算擬合度

$tau

[1] 0.5

$R1

[1] 0.620556

結(jié)果解釋: tau = 0.5的擬合度分別為0.650 556.

第5步: 擬合圖(圖1).

圖1 分位數(shù)回歸的擬合圖. 從上向下的直線依次是tau為0.1, 0.5, 0.9的擬合線.

>plot(income, foodexp, cex = 0.25, type = "n")

>points(income, foodexp, cex = 0.5, col = "blue")

>taus < -c(0.1, 0.5, 0.9)

>for (i in 1:length(taus)) {

>abline(rq(foodexp~income, tau = taus[i]),

lty = 1, col = "green") }

圖2 回歸估計(jì)值的變化圖

第6步: 回歸系數(shù)的圖

>plot(summary(rq(foodexp~income, tau = 5:90/100), se = "boot", bsmethod = "xy"))

圖2是tau = 0.05, 0.51, 0.52…, 0.9的回歸估計(jì)值的變化圖, 圖2(a)是常數(shù)項(xiàng)變化圖, 圖2(b)是系數(shù)變化圖. 水平實(shí)直線表示經(jīng)典回歸的估計(jì)值, 虛直線表示經(jīng)典回歸估計(jì)的95%置信區(qū)間的上下限; 灰色表示分位數(shù)回歸95%置信帶.

6 總結(jié)

相比經(jīng)典的回歸, 分位數(shù)回歸可以描述自變量對于整個(gè)因變量的條件變化影響, 捕捉分布的尾部信息等特點(diǎn). 近年來, 分位數(shù)回歸廣泛應(yīng)用于金融、經(jīng)濟(jì)學(xué)、生物及環(huán)境等領(lǐng)域. 本文主要介紹分位數(shù)回歸的思想、參數(shù)估計(jì)、參數(shù)檢驗(yàn)、模型評價(jià)及其在R中的實(shí)現(xiàn). 最后基于實(shí)例, 詳細(xì)的給出R中實(shí)現(xiàn)的一般步驟.

附錄1 擬似然比檢驗(yàn)R代碼

QLR=function(object, hs=hs, mod=T){

mt<- terms(object)

m <- model.frame(object)

y <- model.response(m)

x < - model.matrix(mt, m, contrasts = object$contrasts)

n < -length(y)

tau<- object$tau

coef<- coefficients(object)

method<-fit$method

p <- dim(x)[2]

rdf < -p - 1

res < -fit$residuals

fit0 < -rq(y ~ 1, method = method, tau = tau)

v1 < -fit$rho

v0 < -fit0$rho

if(mod){

sps < -as.numeric((summary(fit, se = "iid", cov = T, hs = hs))$scale)

sp < -1/sps

}

else{

sps < -as.numeric((summary(fit, se = "ker", cov = T))$scale)

sp < -1/sps

}

f<-2*(v0-v1)/(sp*tau*(1-tau))

Pvalue < -1-pchisq(f, rdf)

structure(list(tau=tau, QLR. Statistic = f, p. value = Pvalue))

}

程序的使用: object是rq返回的對象, hs的參數(shù)設(shè)置與誤差為獨(dú)立同分布中的hs一樣. mod表示計(jì)算稀疏函數(shù)的方法, mod = T使用Koenker改進(jìn)了差分的方法, mod = F使用的核函數(shù)法. 需要注意的是此時(shí)命令只能對應(yīng)一個(gè)tau, 不能同時(shí)對應(yīng)多個(gè)tau.

附錄2 模型擬合度R代碼

R1=function(object){

mt<- terms(object)

m <- model.frame(object)

y <- model.response(m)

tau <- object$tau

coef < - coefficients(object)

method < -fit$method

fit0 < -rq(y ~ 1, method = method, tau = tau)

v1 < -fit$rho

v0 < -fit0$rho

R1 < -1-v1/v0

structure(list(tau = tau, R1 = R1))

}

程序的使用: object是rq返回的對象.

[1] Koenker R, Bassett G. Regression quantiles[J]. Econometrica, 1978, 46(1): 33—50.

[2] 李育安. 分位數(shù)回歸及應(yīng)用簡介[J]. 統(tǒng)計(jì)理論與方法, 2006, 21(3): 35—38.

[3] 荀鵬程, 顧堅(jiān), 顧海燕, 等. 中位數(shù)回歸模型及自回歸模型在北京市SARS發(fā)病預(yù)測中的應(yīng)用[J]. 中國衛(wèi)生統(tǒng)計(jì), 2004, 4: 123—145.

[4] Barrodale I, Roberts F. Solution of anOverdeterminedSystemof Equationsin the l1Norm[J]. Communications of the ACM, 1974, 17(6), 319—320.

[5] 陳建寶, 丁軍軍. 分位數(shù)回歸技術(shù)綜述[J]. 統(tǒng)計(jì)與信息論壇, 2008, 23(3): 89—96.

[6] Portnoy S, Koenker R. The Gaussian Hare and the Laplacian Tortoise: Computability ofSquared-Error Versus Absolute-Error Estimators with Discusssion[J]. Statistical Science, 1997, 12(4): 279—300.

[7] Koenker R. Quantile Regression[M]. London: Cambridge U Press, 2005: 73—81, 202—204.

[8] He X, Hu F. Markov Chain Marginal Bootstrap[J]. Journal of the American Statistical Association, 2002, 97: 783—795.

[9] Parzen M I, Wei L, Ying Z. A resampling method based on pivotal estimating functions[J]. Biometrika, 1994, 81(2), 341—350.

[10] Roger Koenker, Jose A F. Machado. Goodness of Fit and Related Inference Processes for Quantile Regression[J]. Journal of the American Statistical Association, 1999, 94: 1296—1310.

[11] Koenker R, Bassett G. Robust Tests of Heteroscedasticity based on Regression Quantiles[J]. Econometrica, 1982, 50(1), 43—61.

Quantile regressionand itsrealizationin R

HE Feng-xia, WANG Feng-zhu

(School of Mathematics and physics, North China Electric Power University, Beijing 102206, China)

As there isn’t function realizing quasi-likelihood-ratio model test and goodness-of-fit in R, R-codes of quasi-likelihood-ratio and goodness-of-fit programmed are put forward. General process of quantile regression is developed based on an represent-ative example for providing important reference on quantile regression algorithm research.

quantile regression; R; quasi-likelihood-ratio; goodness-of-fit

10.3969/j.issn.1672-6146.2013.03.003

O 212.4

1672-6146(2013)03-0010-04

email: he_fx@126.com.

2013-07-06

(責(zé)任編校: 劉曉霞)