用R语言进行分位数回归.docx

用R语言进行分位数回归.docx

《用R语言进行分位数回归.docx》由会员分享，可在线阅读，更多相关《用R语言进行分位数回归.docx（38页珍藏版）》请在冰点文库上搜索。

用R语言进行分位数回归

用R语言进行分位数回归：

基础篇

詹鹏

（师大学经济管理学院）

本文根据文献资料整理，以介绍方法为主要目的。

作者的主要贡献有：

（1）整理了分位数回归的一些基本原理和方法；

（2）归纳了用R语言处理分位数回归的程序，其中写了两个函数整合估计结果；（3）写了一个分位数分解函数来处理MM2005的分解过程；（4）使用一个数据集进行案例分析，完整地展现了分析过程。

第一节 分位数回归介绍

（一）为什么需要分位数回归？

传统的线性回归模型描述了因变量的条件均值分布受自变量X的影响过程。

其中，最小二乘法是估计回归系数的最基本方法。

如果模型的随机误差项来自均值为零、方差相同的分布，那么回归系数的最小二乘估计为最佳线性无偏估计（BLUE）；如果随机误差项是正态分布，那么回归系数的最小二乘估计与极大似然估计一致，均为最小方差无偏估计（MVUL）。

此时它具有无偏性、有效性等优良性质。

但是在实际的经济生活中，这种假设通常不能够满足。

例如当数据中存在严重的异方差，或后尾、尖峰情况时，最小二乘法的估计将不再具有上述优良性质。

为了弥补普通最小二乘法（OLS）在回归分析中的缺陷，1818年Laplace[2]提出了中位数回归（最小绝对偏差估计）。

在此基础上，1978年Koenker和Bassett[3]把中位数回归推广到了一般的分位数回归（QuantileRegression）上。

分位数回归相对于最小二乘回归，应用条件更加宽松，挖掘的信息更加丰富。

它依据因变量的条件分位数对自变量X进行回归，这样得到了所有分位数下的回归模型。

因此分位数回归相比普通的最小二乘回归，能够更加精确第描述自变量X对因变量Y的变化围，以及条件分布形状的影响。

（二）一个简单的分位数回归模型[4]

假设随机变量的分布函数为

（1）

Y的

分位数的定义为满足

的最小

值，即

（2）

回归分析的基本思想就是使样本值与拟合值之间的距离最短，对于Y的一组随机样本

，样本均值回归是使误差平方和最小，即

                                                                           （3）

样本中位数回归是使误差绝对值之和最小，即

                                                                              （4）

样本分位数回归是使加权误差绝对值之和最小，即

                                          （5）

上式可等价表示为：

其中，

为检查函数（checkfunction），定义为：

其中，

为指示函数（indicatorfunction），z是条件关系式，当z为真时，

；当z为假时，

。

同线性方程y=kx比较，

相当于直线的斜率k，可以看出，

为分段函数，如下图所示。

现假设因变量Y由k个自变量组成的矩阵X线性表示，对于条件均值函数

，通过求解（8）式得到参数估计值

对于条件分位数函数，通过求解（9）式得到参数估计值

式中，

函数表示取函数最小值时

的取值。

（三）分位数回归模型的参数估计算法

1、主要算法

（1）单纯形算法（SimplexMethod）

Koenker和Orey[5]（1993）把分两步解决最优化问题的单纯形算法[6]扩展到所有回归分位数中。

该算法估计出来的参数具有很好的稳定性，但是在处理大型数据时运算的速度会显著降低。

（2）点算法（InteriorPointMethod）

由于单纯形算法在处理大型数据时效率低下，Karmarker提出了点算法[7]；Portnoy和Koenker把这种方法是用在分位数回归中，得出了处理大型数据时点算法的运算速度远快于单纯形算法的结论。

但点算法每计算一步都要进行因数分解，当自变量比较多的时候效率比较低。

其次，如果要达到和单纯形算法一样的精度，就必须进行舍入步骤的计算，者也降低了算法的运行效率。

（3）平滑算法（SmoothingMethod）

上述两种算法都有各自的优点和不足，而有限平滑算法则是一种同时兼顾运算效率以及运算速度的方法。

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:

逻辑参数，是否返回参数估计量的协方差矩阵。

不同参数的结果，可参看下面的程序案例。

（四）分位数分解（MM2005方法）[9]

我们可以进一步运用分位数分解法对各个影响因素进行分解分析[10]。

这里仅介绍MM2005方法。

为讲解方便，这里以各因素对城乡家庭收入的影响为例，观察各个影响因素在不同分位数上对城乡家庭收入差异的影响度的大小。

这里介绍Machado和Mata[11]（2005）提出的分位数分解法，将每个分位数上的城乡收入差异分解为两个部分：

一部分是由于城乡家庭劳动力特征的不同回报率引起的（即分位数回归参数的不同引起的，TheReturnEffects），例如城乡家庭劳动力在相同的教育程度、工作年限以及所处当地的经济发展水平相同的特定因素下不同的回报率引起的家庭人均收入差异；另一部分是由于城乡家庭劳动力的特征变量分布不同引起的（即影响因素变量值的不同引起的，TheCovariateEffect），城乡家庭人均收入这部分的差异会随着样本分布的不同而略有变化。

利用Machado和Mata分位数分解方法的关键是进行反事实分析（thecounter-factualanalysis），我们最关心的一种反事实分析就是，如果城市家庭劳动力按照农村家庭劳动力的分位数回归参数决定家庭人均收入的话，城市家庭的人均收入分布会如何？

这里定义反事实分布为

，其中

表示影响城市家庭人均收入的变量分布，

表示影响农村家庭人均收入的变量在每个分位数上的回归参数。

表示如果城市家庭劳动力按照农村家庭劳动力的分位数回归参数决定家庭人均收入的话，城市家庭的反事实人均收入的大小。

的具体计算步骤为：

（1）确定不同的分位点，分别表示为

。

（2）在农村家庭样本中，分别以

做分位数回归，得到

组分位数回归参数向量

。

（3）将城市家庭样本数据表示为

。

（4）把

（2）中得到的分位数回归参数和（3）中得到得城市家庭子样本变量分布相结合，得到一个新的样本，即反事实分布样本

。

假定在τ分位数下城市家庭人均收入、反事实家庭人均收入和农村家庭人均收入分别为

、

、

。

则不同分位数下的城乡家庭人均收入分布差异可表示为：

等式右边的第一项称为“回报影响（thereturneffect）”，它表示在不同的分位数下，由于城乡家庭劳动力的生产回报率不同所导致的城乡差异部分；等式右边第二项成为“变量影响（thecovariateeffect）”，它表示不同分位数下城乡家庭随机抽样的样本变量分布不同所导致的城乡差异部分。

（五）非线性分位数回归和非参数分位数回归

暂略。

第二节 用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

#1420.1577255.8394

#2541.4117310.9587

#3901.1575485.6800

#4639.0802402.9974

#5750.8756495.5608

#6945.7989633.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

#

#Degreesoffreedom:

235total;233residual

④用summary（）函数可以得到回归模型的详细结果，包括系数和上下限。

#Call:

rq（formula=foodexp~income,tau=0.5,data=engel）

#

#tau:

[1]0.5

#

#Coefficients:

#             coefficientslowerbd upperbd

#（Intercept） 81.48225    53.25915114.01156

#income       0.56018     0.48702  0.60199

⑤coef（）函数得到的系数为向量形式，第一个元素为常数项的系数，第二个及以后为自变量的系数。

⑥summary函数se参数的说明：

A.se=“rank”:

按照Koenker（1994）的排秩方法计算得到的置信区间，默认残差为独立同分布。

注意的是，上下限是不对称的。

Coefficients:

           coefficientslowerbd upperbd

（Intercept） 81.48225    53.25915114.01156

income       0.56018     0.48702  0.60199

B.se=”iid”:

假设残差为独立同分布，用KB（1978）的方法计算得到近似的协方差矩阵。

Coefficients:

           Value   Std.Errortvalue Pr（>|t|）

（Intercept）81.4822513.23908   6.15468 0.00000

income      0.56018 0.01192  46.99766 0.00000

C.se=“nid”:

表示按照Huber方法逼近得到的估计量。

Coefficients:

           Value   Std.Errortvalue Pr（>|t|）

（Intercept）81.4822519.25066   4.23270 0.00003

income      0.56018 0.02828  19.81032 0.00000

D.se=”ker”:

采用Powell（1990）的核估计方法。

Coefficients:

           Value   Std.Errortvalue Pr（>|t|）

（Intercept）81.4822530.21532   2.69672 0.00751

income      0.56018 0.03732  15.01139 0.00000

E.se=”boot”:

采用bootstrap方法自助抽样的方法估计系数的误差标准差。

Coefficients:

           Value   Std.Errortvalue Pr（>|t|）

（Intercept）81.4822525.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.1 99个分位点的系数估计值

图2.2 5个分位点的系数估计值

   2、不同分位点拟合曲线的比较

#散点图

attach（engel）         #打开engel数据集，直接运行其中的列名，就可以调用相应列

plot（income,foodexp,cex=0.25,type="n",  #画图，说明①

    xlab="HouseholdIncome",ylab="FoodExpenditure"）

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（iin1:

length（taus））{      #绘制不同分位点下的拟合直线，颜色为灰色

 abline（rq（foodexp~income,tau=taus[i]）,col="gray"）

}

detach（engel）

图2.3不同分位点下的分位数回归拟合结果比较

   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="FoodExpenditure",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"））

图2.4 10%分位点和90%分位点之间的比较

   上图表示收入（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）

结果：

QuantileRegressionAnalysisofDevianceTable

Model:

foodexp~income

JointTestofEqualityofSlopes:

tauin{ 0.250.50.75 }

 DfResidDfFvalue   Pr（>F）   

1 2     703 15.5572.449e-07***

---

Signif.codes:

 0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

其中P值远小于0.05，故不同分位点下收入对食品支出的影响机制不同。

（五）残差形态的检验

也可以理解为是比较不同分位点的模型之间的关系。

主要有两种模型形式：

（1）位置漂移模型：

不同分位点的估计结果之间的斜率相同或近似，只是截距不同；表现为不同分位点下的拟合曲线是平行的。

（2）位置-尺度漂移模型：

不同分位点的估计结果之间的斜率和截距都不同；表现为不同分位点下的拟合曲线不是平行的。

#残差形态的检验

source（"C:

/ProgramFiles/R/R-2.15.0/library/quantreg/doc/gasprice.R"）

x=gasprice

n=length（x）

p=5

X=cbind（x[（p-1）:

（