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

1、用R语言进行分位数回归用R语言进行分位数回归：基础篇詹鹏（师大学经济管理学院 ）本文根据文献资料整理，以介绍方法为主要目的。作者的主要贡献有：（1）整理了分位数回归的一些基本原理和方法；（2）归纳了用R语言处理分位数回归的程序，其中写了两个函数整合估计结果；（3）写了一个分位数分解函数来处理MM2005的分解过程；（4）使用一个数据集进行案例分析，完整地展现了分析过程。第一节分位数回归介绍（一）为什么需要分位数回归？传统的线性回归模型描述了因变量的条件均值分布受自变量X的影响过程。其中，最小二乘法是估计回归系数的最基本方法。如果模型的随机误差项来自均值为零、方差相同的分布，那么回归系数的最小二

2、乘估计为最佳线性无偏估计（BLUE）；如果随机误差项是正态分布，那么回归系数的最小二乘估计与极大似然估计一致，均为最小方差无偏估计（MVUL）。此时它具有无偏性、有效性等优良性质。但是在实际的经济生活中，这种假设通常不能够满足。例如当数据中存在严重的异方差，或后尾、尖峰情况时，最小二乘法的估计将不再具有上述优良性质。为了弥补普通最小二乘法（OLS）在回归分析中的缺陷，1818年Laplace2提出了中位数回归（最小绝对偏差估计）。在此基础上，1978年Koenker和Bassett3把中位数回归推广到了一般的分位数回归（Quantile Regression）上。分位数回归相对于最小二乘回归，

3、应用条件更加宽松，挖掘的信息更加丰富。它依据因变量的条件分位数对自变量X进行回归，这样得到了所有分位数下的回归模型。因此分位数回归相比普通的最小二乘回归，能够更加精确第描述自变量X对因变量Y的变化围，以及条件分布形状的影响。（二）一个简单的分位数回归模型4假设随机变量的分布函数为 （1）Y的分位数的定义为满足的最小值，即 （2）回归分析的基本思想就是使样本值与拟合值之间的距离最短，对于Y的一组随机样本，样本均值回归是使误差平方和最小，即 （3）样本中位数回归是使误差绝对值之和最小，即 （4）样本分位数回归是使加权误差绝对值之和最小，即 （5）上式可等价表示为： 其中，为检查函数（check f

4、unction），定义为： 其中，为指示函数（indicator function），z是条件关系式，当z为真时，；当z为假时，。同线性方程y=kx比较，相当于直线的斜率k，可以看出，为分段函数，如下图所示。现假设因变量Y由k个自变量组成的矩阵X线性表示，对于条件均值函数，通过求解（8）式得到参数估计值对于条件分位数函数，通过求解（9）式得到参数估计值式中，函数表示取函数最小值时的取值。（三）分位数回归模型的参数估计算法1、主要算法（1）单纯形算法（Simplex Method）Koenker和Orey5（1993）把分两步解决最优化问题的单纯形算法6扩展到所有回归分位数中。该算法估计出来的参

5、数具有很好的稳定性，但是在处理大型数据时运算的速度会显著降低。（2）点算法（Interior Point Method）由于单纯形算法在处理大型数据时效率低下，Karmarker提出了点算法7；Portnoy和Koenker把这种方法是用在分位数回归中，得出了处理大型数据时点算法的运算速度远快于单纯形算法的结论。但点算法每计算一步都要进行因数分解，当自变量比较多的时候效率比较低。其次，如果要达到和单纯形算法一样的精度，就必须进行舍入步骤的计算，者也降低了算法的运行效率。（3）平滑算法（Smoothing Method）上述两种算法都有各自的优点和不足，而有限平滑算法则是一种同时兼顾运算效率以及

6、运算速度的方法。Chen把这种算法扩展到计算回归分位数中8。2、R语言quantreg包中的假设检验加载quantreg包以后，使用summary()函数或summary.rq()函数，可以得到参数系数的一些假设检验统计量。其实，以上两个函数是一致的。在使用summary()的时候，如果sumamry()加载的模型（对象）是分位数回归模型，则会自动调用summary.rq()来处理这个对象。summary.rq()的调用格式为summary(object, se = NULL, covariance=FALSE, hs = TRUE,.)其中主要参数有：# object: 分位数回归对象，根据

7、rq()函数等得到的结果。# se: 用于计算参数估计值标准差的方法，可以选取的值包括：- rank: 根据Koenker(1994)的秩检验得到标准差的估计值。默认情况下假定残差是服从独立同分布。如果补充另一个参数iid=FALSE，则采用Machado(1999)的方法计算标准差（参数的写法：se=”rank”, iid=FALSE）。- iid: （这个与上面提到的iid=FALSE不同，这里是参数se的一个取值，而上面的iid是一个逻辑参数）假定残差服从独立同分布，并按照KB（1978）的方法计算残差。- nid: 用sparsity算法计算的参数估计值标准差。- ker: 用Powe

8、ll(1990)的核密度估计方法得到标准差。- boot: 采用bootstrap自助抽样的方法计算标准差。- 默认情况下，se=NULL且convariance=FALSE，标准差的默认算法是se=”rank”；其他情况下，se默认值为”nid”。# covariance: 逻辑参数，是否返回参数估计量的协方差矩阵。不同参数的结果，可参看下面的程序案例。（四）分位数分解（MM2005方法）9我们可以进一步运用分位数分解法对各个影响因素进行分解分析10。这里仅介绍MM2005方法。为讲解方便，这里以各因素对城乡家庭收入的影响为例，观察各个影响因素在不同分位数上对城乡家庭收入差异的影响度的大小。

9、这里介绍Machado和Mata11（2005）提出的分位数分解法，将每个分位数上的城乡收入差异分解为两个部分：一部分是由于城乡家庭劳动力特征的不同回报率引起的（即分位数回归参数的不同引起的，The Return Effects），例如城乡家庭劳动力在相同的教育程度、工作年限以及所处当地的经济发展水平相同的特定因素下不同的回报率引起的家庭人均收入差异；另一部分是由于城乡家庭劳动力的特征变量分布不同引起的（即影响因素变量值的不同引起的，The Covariate Effect），城乡家庭人均收入这部分的差异会随着样本分布的不同而略有变化。 利用Machado和Mata分位数分解方法的关键是进行反

10、事实分析（the counter-factual analysis），我们最关心的一种反事实分析就是，如果城市家庭劳动力按照农村家庭劳动力的分位数回归参数决定家庭人均收入的话，城市家庭的人均收入分布会如何？这里定义反事实分布为，其中表示影响城市家庭人均收入的变量分布，表示影响农村家庭人均收入的变量在每个分位数上的回归参数。表示如果城市家庭劳动力按照农村家庭劳动力的分位数回归参数决定家庭人均收入的话，城市家庭的反事实人均收入的大小。的具体计算步骤为：（1）确定不同的分位点，分别表示为。（2）在农村家庭样本中，分别以做分位数回归，得到组分位数回归参数向量。 （3）将城市家庭样本数据表示为。（4）把

11、（2）中得到的分位数回归参数和（3）中得到得城市家庭子样本变量分布相结合，得到一个新的样本，即反事实分布样本。假定在分位数下城市家庭人均收入、反事实家庭人均收入和农村家庭人均收入分别为、。则不同分位数下的城乡家庭人均收入分布差异可表示为： 等式右边的第一项称为“回报影响（the return effect）”，它表示在不同的分位数下，由于城乡家庭劳动力的生产回报率不同所导致的城乡差异部分；等式右边第二项成为“变量影响(the covariate effect)”，它表示不同分位数下城乡家庭随机抽样的样本变量分布不同所导致的城乡差异部分。（五）非线性分位数回归和非参数分位数回归暂略。第二节用R语

12、言进行分位数回归（一）安装和加载包 R语言的基本包中没有进行分位数回归的程序包，故需要在官网下载并安装相应的程序包quantreg。在电脑上安装过quantreg包以后，下次不需要再次安装了。但每次使用分位数回归前，需要加载quantreg包。install.package(“quantreg”)# 保持联网的情况下安装包library(“quantreg”) # 加载包help.start() # 进入R帮助首页help(rq) # 获取rq函数的帮助，也可以写成：?rq example(rq) # 显示分位数回归函数rq()的一个简单示例代码（二）一个简单的分位数回归模型及结果data(e

13、ngel) # 加载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）是考察食物支出与家庭收入之

14、间关系的一个数据集，用函数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%分位数下的回归结果。

15、- data=engel指明这里处理的数据集为engel。- method：进行拟合的方法，取值包括：A. 默认值“br”，表示 Barrodale & Roberts 算法的修改版；B. “fn”，针对大数据可以采用的FrischNewton点算法；C. “pfn”，针对特别大数据，使用经过预处理的FrischNewton逼近方法；D. “fnc”，针对被拟合系数特殊的线性不等式约束情况；E. “lasso”和“scad”，基于特定惩罚函数的平滑算法进行拟合。 直接运行fit1，会得到简单的计算结果，如：# Call:# rq(formula = foodexp income, tau =

16、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 bdupper bd # (Intercept)8

17、1.48225 53.25915 114.01156# income 0.56018 0.48702 0.60199 coef()函数得到的系数为向量形式，第一个元素为常数项的系数，第二个及以后为自变量的系数。 summary函数se参数的说明：A. se = “rank”: 按照Koenker(1994)的排秩方法计算得到的置信区间，默认残差为独立同分布。注意的是，上下限是不对称的。Coefficients: coefficients lower bdupper bd (Intercept)81.48225 53.25915 114.01156income 0.56018 0.48702 0

18、.60199B. se=”iid”: 假设残差为独立同分布，用KB（1978）的方法计算得到近似的协方差矩阵。Coefficients: Value Std. Error t valuePr(|t|)(Intercept) 81.48225 13.23908 6.154680.00000income 0.560180.01192 46.997660.00000C. se = “nid”: 表示按照Huber 方法逼近得到的估计量。Coefficients: Value Std. Error t valuePr(|t|)(Intercept) 81.48225 19.25066 4.232700

19、.00003income 0.560180.02828 19.810320.00000D. se=”ker”: 采用Powell(1990)的核估计方法。Coefficients: Value Std. Error t valuePr(|t|)(Intercept) 81.48225 30.21532 2.696720.00751income 0.560180.03732 15.011390.00000E. se=”boot”: 采用bootstrap方法自助抽样的方法估计系数的误差标准差。Coefficients: Value Std. Error t valuePr(|t|)(Interc

20、ept) 81.48225 25.23647 3.228750.00142income 0.560180.03194 17.537520.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) # 新

21、建一个图形窗口，可以去掉这句plot(fit2)结果：图2.199个分位点的系数估计值图2.25个分位点的系数估计值 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.5abline( rq(foodexp income, tau=0.5

22、), 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=tausi), col=gray )detach(engel)图2.3 不同分位点下的分位数回归拟合结果比较 3、穷人和富人的消费分布比较# 比较穷人（收入在10%分位点的那个

23、人）和富人（收入在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$sol1, # 每个分位点的tau值qs.poor = c( c(1,x.poor) %*% z$sol4:5, )# 10%分位点的收入的消费估计值qs.rich = c( c(

24、1,x.rich) %*% z$sol4: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.poor1,qs.poor), do.points=F, add=T)plot(stepfun(ps,c(qs.poor1,qs.rich), do.points=F, add=T

25、, col.hor=gray, col.vert=gray)ps.wts = ( c(0,diff(ps) + c(diff(ps),0) )/2ap = 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=blac

26、k)legend(topright, c(poor,rich), lty=c(1,1), col=c(black,gray)图2.410%分位点和90%分位点之间的比较 上图表示收入（income）为10%分位点处（poor，穷人）和90%分位点处（rich，富人）的食品支出的比较。从左图可以发现，对于穷人而言，在不同分位点估计的食品消费差别不大。而对于富人而言，在不同分位点对食品消费的差别比较大。右图反应了穷人和富人的食品消费分布曲线。穷人的食品消费集中于400左右，比较陡峭；而富人的消费支出集中于800到1200之间，比较分散。（四）模型比较# 比较不同分位点下，收入对食品支出的影响机制是

27、否相同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 TableModel: foodexp incomeJoint Test of Equality of Slopes: tau in 0.25 0.5 0.75Df Resid Df F value Pr(F) 12 70315.557 2.4

28、49e-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 = gaspricen = length(x)p = 5X = cbind(x(p-1):(