本文根据文献资料整理，以介绍方法为主要目的。作者的主要贡献有：（1）整理了分位数回归的一些基本原理和方法；（2）归纳了用R语言处理分位数回归的程序，其中写了两个函数整合估计结果；（3）使用一个数据集进行案例分析，完整地展现了分析过程。
第一节 分位数回归介绍
（一）为什么需要分位数回归？
传统的线性回归模型描述了因变量的条件均值分布受自变量X的影响过程。其中，最小二乘法是估计回归系数的最基本方法。如果模型的随机误差项来自均值为零、方差相同的分布，那么回归系数的最小二乘估计为最佳线性无偏估计（BLUE）；如果随机误差项是正态分布，那么回归系数的最小二乘估计与极大似然估计一致，均为最小方差无偏估计（MVUL）。此时它具有无偏性、有效性等优良性质。
但是在实际的经济生活中，这种假设通常不能够满足。例如当数据中存在严重的异方差，或后尾、尖峰情况时，最小二乘法的估计将不再具有上述优良性质。为了弥补普通最小二乘法（OLS）在回归分析中的缺陷，1818年Laplace[2]提出了中位数回归（最小绝对偏差估计）。在此基础上，1978年Koenker和Bassett[3]把中位数回归推广到了一般的分位数回归（Quantile Regression）上。
分位数回归相对于最小二乘回归，应用条件更加宽松，挖掘的信息更加丰富。它依据因变量的条件分位数对自变量X进行回归，这样得到了所有分位数下的回归模型。因此分位数回归相比普通的最小二乘回归，能够更加精确第描述自变量X对因变量Y的变化范围，以及条件分布形状的影响。
（二）一个简单的分位数回归模型[4]
假设随机变量的分布函数为
Y的分位数的定义为满足
回归分析的基本思想就是使样本值与拟合值之间的距离最短，对于Y的一组随机样本
（3）
样本中位数回归是使误差绝对值之和最小
（4）
样本分位数回归是使加权误差绝对值之和最小，即
上式可等价表示为：
其中，为检查函数（check function），定义为：
其中，
现假设因变量Y由k个自变量组成的矩阵X线性表示，对于条件均值函数
对于条件分位数函数，通过求解（7）式得到参数估计值
式中，
（三）分位数回归模型的参数估计算法
1、主要算法
（1）单纯形算法（Simplex Method）
Koenker和Orey[5]（1993）把分两步解决最优化问题的单纯形算法[6]扩展到所有回归分位数中。该算法估计出来的参数具有很好的稳定性，但是在处理大型数据时运算的速度会显著降低。
（2）内点算法（Interior Point Method）
由于单纯形算法在处理大型数据时效率低下，Karmarker提出了内点算法[7]；Portnoy和Koenker把这种方法是用在分位数回归中，得出了处理大型数据时内点算法的运算速度远快于单纯形算法的结论。但内点算法每计算一步都要进行因数分解，当自变量比较多的时候效率比较低。其次，如果要达到和单纯形算法一样的精度，就必须进行舍入步骤的计算，者也降低了算法的运行效率。
（3）平滑算法（Smoothing Method）
上述两种算法都有各自的优点和不足，而有限平滑算法则是一种同时兼顾运算效率以及运算速度的方法。Chen把这种算法扩展到计算回归分位数中[8]。
2、R语言quantreg包中的假设检验
加载quantreg包以后，使用summary()函数或summary.rq()函数，可以得到参数系数的一些假设检验统计量。其实，以上两个函数是一致的。在使用summary()的时候，如果sumamry()加载的模型（对象）是分位数回归模型，则会自动调用summary.rq()来处理这个对象。summary.rq()的调用格式为
summary(object, se = NULL, covariance=FALSE, hs = TRUE, ...)
其中主要参数有：
# object: 分位数回归对象，根据rq()函数等得到的结果。
# se: 用于计算参数估计值标准差的方法，可以选取的值包括：
- rank: 根据Koenker(1994)的秩检验得到标准差的估计值。默认情况下假定残差是服从独立同分布。如果补充另一个参数iid=FALSE，则采用Machado(1999)的方法计算标准差（参数的写法：se=”rank”, iid=FALSE）。
- iid: （这个与上面提到的iid=FALSE不同，这里是参数se的一个取值，而上面的iid是一个逻辑参数）假定残差服从独立同分布，并按照KB（1978）的方法计算残差。
- nid: 用sparsity算法计算的参数估计值标准差。
- ker: 用Powell(1990)的核密度估计方法得到标准差。
- boot: 采用bootstrap自助抽样的方法计算标准差。
- 默认情况下，se=NULL且convariance=FALSE，标准差的默认算法是se=”rank”；其他情况下，se默认值为”nid”。
# covariance: 逻辑参数，是否返回参数估计量的协方差矩阵。
不同参数的结果，可参看下面的程序案例。
第二节 用R语言进行分位数回归
（一）安装和加载包
R语言的基本包中没有进行分位数回归的程序包，故需要在官网下载并安装相应的程序包quantreg。在电脑上安装过quantreg包以后，下次不需要再次安装了。但每次使用分位数回归前，需要加载quantreg包。
install.package(“quantreg”) # 保持联网的情况下安装包
library(“quantreg”) # 加载包
help.start() # 进入R帮助首页
help(rq) # 获取rq函数的帮助，也可以写成：?rq
example(rq) # 显示分位数回归函数rq()的一个简单示例代码
（二）一个简单的分位数回归模型及结果
data(engel) # 加载quantreg包自带的数据集，见说明①
fit1 = rq(foodexp ~ income, tau = 0.5, data = engel) # 进行分位数回归，见说明②
fit1 # 直接显示分位数回归的模型和系数，见说明③
summary(fit1) # 得到更加详细的显示结果，见说明④
r1 = resid(fit1) # 得到残差序列，并赋值为变量r1
c1 = coef(fit1) # 得到模型的系数，并赋值给变量c1，见说明⑤
summary(fit1, se = “nid”) # 通过设置参数se，可以得到系数的假设检验，说明⑥
说明：① engel（1857）是考察食物支出与家庭收入之间关系的一个数据集，用函数head(engel)可以查看前六行的值：
# income foodexp
# 1 420.1577 255.8394
# 2 541.4117 310.9587
# 3 901.1575 485.6800
# 4 639.0802 402.9974
# 5 750.8756 495.5608
# 6 945.7989 633.7978
② 这里因变量为foodexp，即食物支出。自变量为income，即家庭收入。
- tau表示计算50%分位点的参数，这里可以同时计算多个分位点的分位数回归结果，如tau=c(0.1,0.5,0.9)是同时计算10%、50%、90%分位数下的回归结果。
- data=engel指明这里处理的数据集为engel。
- method：进行拟合的方法，取值包括：A. 默认值“br”，表示 Barrodale & Roberts 算法的修改版；B. “fn”，针对大数据可以采用的Frisch–Newton内点算法；C. “pfn”，针对特别大数据，使用经过预处理的Frisch–Newton逼近方法；D. “fnc”，针对被拟合系数特殊的线性不等式约束情况；E. “lasso”和“scad”，基于特定惩罚函数的平滑算法进行拟合。
③ 直接运行fit1，会得到简单的计算结果，如：
# Call:
# rq(formula = foodexp ~ income, tau = 0.5, data = engel)
#
# Coefficients:
# (Intercept) income
# 81.4822474 0.5601806
#
# Degrees of freedom: 235 total; 233 residual
④ 用summary()函数可以得到回归模型的详细结果，包括系数和上下限。
# Call: rq(formula = foodexp ~ income, tau = 0.5, data = engel)
#
# tau: [1] 0.5
#
# Coefficients:
# coefficients lower bd upper bd
# (Intercept) 81.48225 53.25915 114.01156
# income 0.56018 0.48702 0.60199
⑤ coef()函数得到的系数为向量形式，第一个元素为常数项的系数，第二个及以后为自变量的系数。
⑥ summary函数se参数的说明：
A. se = “rank”: 按照Koenker(1994)的排秩方法计算得到的置信区间，默认残差为独立同分布。注意的是，上下限是不对称的。
Coefficients:
coefficients lower bd upper bd
(Intercept) 81.48225 53.25915 114.01156
income 0.56018 0.48702 0.60199
B. se=”iid”: 假设残差为独立同分布，用KB（1978）的方法计算得到近似的协方差矩阵。
Coefficients:
Value Std. Error t value Pr(>|t|)
(Intercept) 81.48225 13.23908 6.15468 0.00000
income 0.56018 0.01192 46.99766 0.00000
C. se = “nid”: 表示按照Huber 方法逼近得到的估计量。
Coefficients:
Value Std. Error t value Pr(>|t|)
(Intercept) 81.48225 19.25066 4.23270 0.00003
income 0.56018 0.02828 19.81032 0.00000
D. se=”ker”: 采用Powell(1990)的核估计方法。
Coefficients:
Value Std. Error t value Pr(>|t|)
(Intercept) 81.48225 30.21532 2.69672 0.00751
income 0.56018 0.03732 15.01139 0.00000
E. se=”boot”: 采用bootstrap方法自助抽样的方法估计系数的误差标准差。
Coefficients:
Value Std. Error t value Pr(>|t|)
(Intercept) 81.48225 25.23647 3.22875 0.00142
income 0.56018 0.03194 17.53752 0.00000
（三）不同分位点下的回归结果比较
1、不同分为点系数估计值的比较
# 不同分位点下的系数估计值的比较
fit1 = summary( rq(foodexp ~ income, tau = 2:98/100) )
fit2 = summary( rq(foodexp ~ income, tau = c(0.05,0.25,0.5,0.75,0.95)) )
windows(5,5) # 新建一个图形窗口，可以去掉这句
plot(fit1)
windows(5,5) # 新建一个图形窗口，可以去掉这句
plot(fit2)
结果：
2、不同分位点拟合曲线的比较
# 散点图
attach(engel) # 打开engel数据集，直接运行其中的列名，就可以调用相应列
plot(income,foodexp,cex=0.25,type='n', # 画图，说明①
xlab='Household Income', ylab='Food Expenditure')
points(income,foodexp,cex=0.5,col='blue') # 添加点，点的大小为0.5
abline( rq(foodexp ~ income, tau=0.5), col='blue' ) # 画中位数回归的拟合直线，颜色蓝
abline( lm(foodexp ~ income), lty = 2, col='red' ) # 画普通最小二乘法拟合直线，颜色红
taus = c(0.05, 0.1, 0.25, 0.75, 0.9, 0.95)
for(i in 1:length(taus)){ # 绘制不同分位点下的拟合直线，颜色为灰色
abline( rq(foodexp ~ income, tau=taus[i]), col='gray' )
}
detach(engel)
3、穷人和富人的消费分布比较
# 比较穷人（收入在10%分位点的那个人）和富人（收入在90%分位点的那个人）的估计结果
# rq函数中，tau不在[0,1]时，表示按最细的分位点划分方式得到分位点序列
z = rq(foodexp ~ income, tau=-1)
z$sol # 这里包含了每个分位点下的系数估计结果
x.poor = quantile(income, 0.1) # 10%分位点的收入
x.rich = quantile(income, 0.9) # 90%分位点的收入
ps = z$sol[1,] # 每个分位点的tau值
qs.poor = c( c(1,x.poor) %*% z$sol[4:5,] ) # 10%分位点的收入的消费估计值
qs.rich = c( c(1,x.rich) %*% z$sol[4:5,] ) # 90%分位点的收入的消费估计值
windows(, 10,5)
par(mfrow=c(1,2)) # 把绘图区域划分为一行两列
plot(c(ps,ps),c(qs.poor,qs.rich),type='n', # type=”n”表示初始化图形区域，但不画图
xlab=expression(tau), ylab='quantile')
plot(stepfun(ps,c(qs.poor[1],qs.poor)), do.points=F,
add=T)
plot(stepfun(ps,c(qs.poor[1],qs.rich)), do.points=F,
add=T, col.hor='gray', col.vert='gray')
ps.wts = ( c(0,diff(ps)) + c(diff(ps),0) )/2
ap = akj(qs.poor, z=qs.poor, p=ps.wts)
ar = akj(qs.rich, z=qs.rich, p=ps.wts)
plot(c(qs.poor,qs.rich), c(ap$dens, ar$dens),
type='n', xlab='Food Expenditure', ylab='Density')
lines(qs.rich,ar$dens,col='gray')
lines(qs.poor,ap$dens,col='black')
legend('topright', c('poor','rich'), lty=c(1,1),
col=c('black','gray'))
上图表示收入（income）为10%分位点处（poor，穷人）和90%分位点处（rich，富人）的食品支出的比较。从左图可以发现，对于穷人而言，在不同分位点估计的食品消费差别不大。而对于富人而言，在不同分位点对食品消费的差别比较大。右图反应了穷人和富人的食品消费分布曲线。穷人的食品消费集中于400左右，比较陡峭；而富人的消费支出集中于800到1200之间，比较分散。
（四）模型比较
# 比较不同分位点下，收入对食品支出的影响机制是否相同
fit1 = rq(foodexp ~ income, tau = 0.25)
fit2 = rq(foodexp ~ income, tau = 0.5)
fit3 = rq(foodexp ~ income, tau = 0.75)
anova(fit1,fit2,fit3)
结果：
Quantile Regression Analysis of Deviance Table
Model: foodexp ~ income
Joint Test of Equality of Slopes: tau in { 0.25 0.5 0.75 }
Df Resid Df F value Pr(>F)
1 2 703 15.557 2.449e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
其中P值远小于0.05，故不同分位点下收入对食品支出的影响机制不同。
（五）残差形态的检验
也可以理解为是比较不同分位点的模型之间的关系。主要有两种模型形式：
（1）位置漂移模型：不同分位点的估计结果之间的斜率相同或近似，只是截距不同；表现为不同分位点下的拟合曲线是平行的。
（2）位置-尺度漂移模型：不同分位点的估计结果之间的斜率和截距都不同；表现为不同分位点下的拟合曲线不是平行的。
# 残差形态的检验
source('C:/Program Files/R/R-2.15.0/library/quantreg/doc/gasprice.R')
x = gasprice
n = length(x)
p = 5
X = cbind(x[(p-1):(n-1)],x[(p-2):(n-2)],x[(p-3):(n-3)],x[(p-4):(n-4)])
y = x[p:n]
# 位置漂移模型的检验
T1 = KhmaladzeTest(y~X, taus = -1, nullH='location')
T2 = KhmaladzeTest(y~X, taus = 10:290/300,
nullH='location', se='ker')
# 位置尺度漂移模型的检验
T3 = KhmaladzeTest(y~X, taus = -1, nullH='location-scale')
T4 = KhmaladzeTest(y~X, taus = 10:290/300,
nullH='location-scale', se='ker')
结果：运行T1，可以查看其检验结果。其中nullH表示原假设为“location”，即原假设为位置漂移模型。Tn表示模型整体的检验，统计量为4.8。THn是对每个自变量的检验。比较T1和T3的结果（T3的原假设为“位置尺度漂移模型”），T1的统计量大于T3的统计量，可见相对而言，拒绝“位置漂移模型”的概率更大，故相对而言“位置尺度漂移模型”更加合适一些。
> T1
$nullH
[1] 'location'
$Tn
[1] 4.803762
$THn
X1 X2 X3 X4
1.0003199 0.5321693 0.5020834 0.8926828
attr(,'class')
[1] 'KhmaladzeTest'
> T3
$nullH
[1] 'location-scale'
$Tn
[1] 2.705583
$THn
X1 X2 X3 X4
1.2102899 0.6931785 0.5045163 0.8957127
attr(,'class')
[1] 'KhmaladzeTest'
（六）非线性分位数回归
这里的非线性函数为Frank copula函数。
## Demo of nonlinear quantile regression model based on Frank copula
vFrank <- function(x,="" df,="" delta,="" u)="" #="">某个非线性过程，得到的是[0,1]的值
-log(1-(1-exp(-delta))/(1+exp(-delta*pt(x,df))*((1/u)-1)))/delta
# 非线性模型
FrankModel <- function(x,="" delta,="" mu,sigma,="" df,="" tau)="" {="">
z <- qt(vfrank(x,="" df,="" delta,="" u="tau),">
mu + sigma*z
}
n <- 200="" #="">样本量
df <- 8="" #="">自由度
delta <- 8="" #="">初始参数
set.seed(1989)
x <- sort(rt(n,df))="" #="">生成基于T分布的随机数
v <- vfrank(x,="" df,="" delta,="" u="runif(n))" #="">基于x生成理论上的非参数对应值
y <- qt(v,="" df)="" #="" v="">对应的T分布统计量
windows(5,5)
plot(x, y, pch='o', col='blue', cex = .25) # 散点图
Dat <- data.frame(x="x," y="y)" #="">基本数据集
us <->
for(i in 1:length(us)){
v <- vfrank(x,="" df,="" delta,="" u="">
lines(x, qt(v,df)) # v为概率，计算每个概率对应的T分布统计量
}
cfMat <- matrix(0,="" 3,="" length(us)+1)="" #="">初始矩阵，用于保存结果的系数
for(i in 1:length(us)) {
tau <->
cat('tau = ', format(tau), '.. ')
fit <->nlrq(y ~ FrankModel(x, delta,mu,sigma, df = 8, tau = tau), # 非参数模型
data = Dat, tau = tau, # data表明数据集,tau分位数回归的分位点
start= list(delta=5, mu = 0, sigma = 1), # 初始值
trace = T) # 每次运行后是否把结果显示出来
lines(x, predict(fit, newdata=x), lty=2, col='red') # 绘制预测曲线
cfMat[i,1] <- tau="" #="">保存分位点的值
cfMat[i,2:4] <- coef(fit)="" #="">保存系数到cfMat矩阵的第i行
cat('\n') # 如果前面把每步的结果显示出来，则每次的结果之间添加换行符
}
colnames(cfMat) <->分位点',names(coef(fit))) # 给保存系数的矩阵添加列名
cfMat
结果：
拟合结果：（过程略）
> cfMat
分位点 delta mu sigma
[1,] 0.25 14.87165 -0.20530041 0.9134657
[2,] 0.50 16.25362 0.03232525 0.9638209
[3,] 0.75 12.09836 0.11998614 0.9423476
（七）半参数和非参数分位数回归
非参数分位数回归在局部多项式的框架下操作起来更加方便。可以基于以下函数。
# 2-7-1 半参数模型----
fit<>
# 其中bs()表示按b-spline的非参数拟合
# 2-7-2 非参数方法
lprq <>
# 这是自定义的一个非参数计算函数，在其他数据下同样可以使用
xx<-seq(min(x),max(x),length=m) #="">个监测点
fv<>
dv<>
for(i in 1:length(xx)){
z<-x-xx[i]>
wx<-dnorm(z )="" #="">核函数为正态分布,dnorm计算标准正态分布的密度值
r<-rq(y~z,weights=wx,tau=tau, #="">上面计算得到的密度值为权重
ci=FALSE)
fv[i]<>
dv[i]<>
}
list(xx=xx,fv=fv,dv=dv) # 输出结果
}
library(MASS)
data(mcycle)
attach(mcycle)
windows(5,5) # 非参数的结果一般是通过画图查看的
plot(times,accel,xlab='milliseconds',ylab='acceleration')
hs<-c(1,2,3,4) #="">选择不同窗宽进行估计
for(i in hs){
h=hs[i]
fit<-lprq(times,accel,h=h,tau=0.5) #="">关键拟合函数
lines(fit$xx,fit$fv,lty=i)
}
legend(45,-70,c('h=1','h=2','h=3','h=4'),
lty=1:length(hs))
# 2-7-3 非参数回归的另一个方法----
# 考察最大的跑步速度与体重的关系
data(Mammals)
attach(Mammals)
x<-log(weight) #="">取得自变量的值
y<-log(speed) #="">取得因变量的值
plot(x,y,xlab='Weightinlog(Kg)',ylab='Speedinlog(Km/hour)',
type='n')
points(x[hoppers],y[hoppers],pch='h',col='red')
points(x[specials],y[specials],pch='s',col='blue')
others<>
points(x[others],y[others],col='black',cex=0.75)
fit<-rqss(y~qss(x,lambda=1),tau=0.9) #="">关键的拟合步骤
plot(fit,add=T) # add=T表示在上图中添加这里的曲线
（八）分位数回归的分解
# MM2005分位数分解的函数，改变参数可直接使用
MM2005 = function(formu,taus, data, group, pic=F){
# furmu 为方程，如foodexp~income
# taus 为不同的分位数
# data 总的数据集
# group 分组指标，是一个向量，用于按行区分data
# pic 是否画图，如果分位数比较多，建议不画图
engel1 = data[group==1,]
engel2 = data[group==2,]
# 开始进行分解
fita = summary( rq(formu, tau = taus, data = engel1 ) )
fitb = summary( rq(formu, tau = taus, data = engel2 ) )
tab = matrix(0,length(taus),4)
colnames(tab) = c('分位数','总差异','回报影响','变量影响')
rownames(tab) = rep('',dim(tab)[1])
for( i in 1:length(taus) ){
ya = cbind(1,engel1[,names(engel1)!=formu[[2]]] ) %*% fita[[i]]$coef[,1]
yb = cbind(1,engel2[,names(engel2)!=formu[[2]]] ) %*% fitb[[i]]$coef[,1]
# 这里以group==1为基准模型，用group==2的数据计算反常规模型拟合值
ystar = cbind(1,engel2[,names(engel2)!=formu[[2]]] ) %*% fita[[i]]$coef[,1]
ya = mean(ya)
yb = mean(yb)
ystar = mean(ystar)
tab[i,1] = fita[[i]]$tau
tab[i,2] = yb - ya
tab[i,3] = yb - ystar # 回报影响，数据相同，模型不同：模型机制的不同所产生的差异
tab[i,4] = ystar - ya # 变量影响，数据不同，模型相同：样本点不同产生的差异
}
# 画图
if( pic ){
attach(engel)
windows(5,5)
plot(income, foodexp, cex=0.5, type='n',main='两组分位数回归结果比较')
points(engel1, cex=0.5, col=2)
points(engel2, cex=0.5, col=3)
for( i in 1:length(taus) ){
abline( fita[[i]], col=2 )
abline( fitb[[i]], col=3 )
}
detach(engel)
}
# 输出结果
tab
}
# 下面是用一个样本数据做的测试
data(engel)
group = c(rep(1,100),rep(2,135)) # 取前100个为第一组，后135个第二组
taus = c(0.05,0.25,0.5,0.75,0.95) # 需要考察的不同分位点
MM2005(foodexp~income, taus, data = engel, group=group, pic=F) # 参数说明，见①
说明：①自编函数MM2005的参数说明：
函数调用形式MM2005 (formu, taus, data, group, pic=F)，其中
# furmu 为回归方程，如foodexp~income；
# taus 为不同的分位数，如taus=c(0.05,0.5,0.95)；
# data 总的数据集，如上例中的engel数据集；
# group 分组指标，是一个向量，用于按行区分data，第一组为1，第二组为2；目前仅能分两组；
# pic 逻辑参数：是否画图。如果分位数比较多，建议不画图；默认不画图，pic=F；如果想画图，则添加参数pic=T。
②最终结果：
> MM2005(foodexp~income, taus, data = engel, group=group, pic=F)
分位数 总差异 回报影响 变量影响
0.05 -30.452061 -72.35939 41.90733
0.25 -2.017317 -46.20125 44.18394
0.50 30.941212 -23.24042 54.18163
0.75 43.729025 -15.76283 59.49185
0.95 52.778985 -11.29932 64.07830
联系客服
