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

何鳳霞, 王鳳竹

?

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

何鳳霞*, 王鳳竹

(華北電力大學 數理學院, 北京, 102206)

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

分位數回歸; R; 擬似然比; 擬合度

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

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

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

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

1 分位數回歸的簡介

2 參數的估計

2.1 單純形算法

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

2.2 內點法

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

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

3 參數檢驗

本文主要介紹參數檢驗常用方法的3種: 誤差為獨立同分布、誤差為非獨立同分布、Bootstrap法.

3.1 誤差為獨立同分布

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

3.2 誤差為非獨立同分布

對于誤差為非獨立同分布時, 參數的極限分布比較復雜, 如位移尺度假設模型:

3.3 Bootstrap法

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

3.3.1-成對Bootstrap

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

3.32 馬爾可夫鏈邊際Bootstrap

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

4 模型檢驗及擬合度

4.1 模型的檢驗

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

4.2 模型的擬合度

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

5 實例

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

第1步: 參數估計.

采用單純形法求解參數估計tau = 0.5的回歸系數. 命令:

>library(quantreg)#調用quantreg包

>data(engel)#讀取engel數據

>attach(engel)#連接數據框的變量名到內存中, 以便直接調用變量

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

>fit#查看回歸估計參數

Coefficients:

Intercept income

81.482 0.560

結果解釋: 中位數回歸的常數項為81.482, 系數為0.56.

第2步: 參數的檢驗.

采用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

結果解釋: 可以得到常數項的= 0.003 < 0.05, 系數的= 0.000 < 0.05, 所以在顯著水平為0.05下參數均顯著;

第3步: 模型的檢驗.

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

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

$tau

[1] 0.5

$QLR.statistic

[1] 437.4879

$p.value

[1] 0

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

第4步: 擬合度的計算.

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

>R1(fit)#計算擬合度

$tau

[1] 0.5

$R1

[1] 0.620556

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

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

圖1 分位數回歸的擬合圖. 從上向下的直線依次是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 回歸估計值的變化圖

第6步: 回歸系數的圖

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

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

6 總結

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

附錄1 擬似然比檢驗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的參數設置與誤差為獨立同分布中的hs一樣. mod表示計算稀疏函數的方法, mod = T使用Koenker改進了差分的方法, mod = F使用的核函數法. 需要注意的是此時命令只能對應一個tau, 不能同時對應多個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] 李育安. 分位數回歸及應用簡介[J]. 統(tǒng)計理論與方法, 2006, 21(3): 35—38.

[3] 荀鵬程, 顧堅, 顧海燕, 等. 中位數回歸模型及自回歸模型在北京市SARS發(fā)病預測中的應用[J]. 中國衛(wèi)生統(tǒng)計, 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] 陳建寶, 丁軍軍. 分位數回歸技術綜述[J]. 統(tǒng)計與信息論壇, 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

(責任編校: 劉曉霞)