用R语言进行分位数回归

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

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

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

4、#39;(n)为指示函数(indicator function ), z是条件关系式，当z为真时，'(")1 ；当z为假时，'(之)口 o同线性方程y二kx比为分段函数，如下图所示。较，一(* 町相当于直线的斜率k,可以看出，80图】检宣函数示意图现假设因变量Y由k个自变量组成的矩阵X线性表示，对于条件均值函数国门户通过求解式得到参数估计值方=argmiiix (2 一切。对于条件分位数函数，通过求解(9)式得到参数估计值& =£ 而！2 (工 NCk 一格力)j-i式中， '曜ndn-函数表示取函数最小值时的取值。(三)分位数回归模型的参数

5、估计算法1、主要算法(1)单纯形算法(Simplex Method )Koenker和Orey5( 1993)把分两步解决最优化问题的单纯形算法 6扩展到所有回归分位数中。该算法估计出来的参数具有很好的稳定性，但是在处理大型数据时运算的速度会显著降低。(2)内点算法(Interior Point Method )由于单纯形算法在处理大型数据时效率低下，Karmarker提出了内点算法7 ; Portnoy和Koenker把这种方法是用在分位数回归中，得出了处理大型 数据时内点算法的运算速度远快于单纯形算法的结论。但内点算法每计算一步都要进行因数分解，当自变量比较多的时候效率比较低。其次，如果要

6、达到和单 纯形算法一样的精度，就必须进行舍入步骤的计算，者也降低了算法的运行效率。(3)平滑算法(Smoothing Method)Chen把这种算法扩展到计算回归分位数中上述两种算法都有各自的优点和不足，而有限平滑算法则是一种同时兼顾运算效率以及运算速度的方法。82、R语言quantreg包中的假设检验加载quantreg包以后，使用summary。函数或summary.rq()函数，可以得到参数系数的一些假设检验统计量。其实，以上两个函数是一致的。在使用summary。的时候，如果sumamry()加载的模型(对象)是分位数回归模型，则会自动调用 summary.rq()来处理这个对象。s

7、ummary.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不同，这

8、里是参数se的一个取值，而上面的iid是一个逻辑参数)假定残差服从独立同分布，并按照KB 1978) 的方法计算残差。nid:用sparsity 算法计算的参数估计值标准差。ker:用Powell(1990)的核密度估计方法得到标准差。boot:采用bootstrap自助抽样的方法计算标准差。默认情况下，se=NULLS, convariance=FALSE,标准差的默认算法是se二“ rank” ；其他情况下，se默认值为“ nid ”。# covariance: 逻辑参数，是否返回参数估计量的协方差矩阵。 不同参数的结果，可参看下面的程序案例。(四)分位数分解(MM2005?法)9我们可以

9、进一步运用分位数分解法对各个影响因素进行分解分析10。这里仅介绍MM2005?法。为讲解方便，这里以各因素对城乡家庭收入的影响为例，观察各个影响因素在不同分位数上对城乡家庭收入差异的影响度的大小。这里介绍Machado?口Mata11 (2005)提出的分位数分解法，将每个分位数上的城乡收入差异分解为两个部分：一部分是由于城乡家庭劳动力特征的不同回报率引起的(即分位数回归参数的不同引起的，The Return Effects ),例如城乡家庭劳动力在相同的教育程度、工作年限以及所处当地的经济发展水平相同的特定因素下不同的 回报率引起的家庭人均收入差异；另一部分是由于城乡家庭劳动力的特征变量分布

10、不同引起的(即影响因素变量值的不同引起的，The Covariate Effect ),城乡家庭人均收入这部分的差异会随着样本分布的不同而略有变化。利用Machado和Mata分位数分解方法的关键是进行反事实分析(the counter-factual analysis),我们最关心的一种反事实分析就是，如果 城市家 庭劳动力按照农村家庭劳动力的分位数回归参数决定家庭人均收入的话，城市家庭的人均收入分布会如何？这里定义反事实分布为F" |N1P其中不表示影响城市家庭人均收入的变量分布,表示影响农村家庭人均收入的变量在每个分位数上的回归参数。表示如果城市家庭劳动力按的具体计算步骤为：(

11、1)确定不同的分位点，分别表示为0 (2)在农村家庭样本中，分别以 FP-做分位数回归，得到组分位数回归参数向量。(3)将城市照农村家庭劳动力的分位数回归参数决定家庭人均收入的话，城市家庭的反事实人均收入的大小。家庭样本数据表示为(4)把(2)中得到的分位数回归参数和(3)中得到得城市家庭子样本变量分布相结合，得到一个新的样本，即反事实分布样。则不同分位数下的城乡家庭人均收入假定在T分位数下城市家庭人均收入、反事实家庭人均收入和农村家庭人均收入分别为分布差异可表示为：J - M 二 Of>等式右边的第一项称为"回报影响(the return effect )”,它表示在不同的分

12、位数下，由于城乡家庭劳动力的生产回报率不同所导致的城乡差异部分；等式右边第二项成为"变量影响(the covariate effect) ”,它表示不同分位数下城乡家庭随机抽样的样本变量分布不同所导致的城乡差异部分。(五)非线性分位数回归和非参数分位数回归暂略。第二节 用R语言进行分位数回归(一)安装和加载包R语言的基本包中没有进行分位数回归的程序包，故需要在官网下载并安装相应的程序包quantreg。在电脑上安装过quantreg包以后，下次不需要再次安装了。但每次使用分位数回归前，需要加载quantreg包。install.package(quantreg ")#保持联

13、网的情况下安装包library( "quantreg " )加载包help.start()进入R帮助首页help(rq)获取rq函数的帮助，也可以写成：？rqexample(rq)显示分位数回归函数rq()的一个简单示例代码(二)一个简单的分位数回归模型及结果data(engel)加载quantreg包自带的数据集，见说明 fitl = rq(foodexp income, tau = 0.5, data = engel) #进行分位数回归,见说明fit1直接显示分位数回归的模型和系数，见说明 summary(fit1)得到更加详细的显示结果，见说明r1 = resid(f

14、it1)得到残差序列，并赋值为变量r1c1 = 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

15、633.7978 这里因变量为foodexp,即食物支出。自变量为income,即家庭收入。- tau表示af算50%H立点的参数，这里可以同时计算多个分位点的分位数回归结果，如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",针对特别大数据，使用经过预处

16、理的 Frisch Newton 逼近方法；D. “fnc”，针对被拟合系数特殊的线性不等式约束情况；E. "lasso”和“scad”，基于特定惩罚函数的平滑算法进行拟合。 直接运行力廿，会得到简单的计算结果，如： # Call: # rq(formula = foodexp income, tau = 0.5, data = engel) # # Coefficients: # (Intercept) income # 81.48224740.5601806# # Degrees of freedom: 235 total; 233 residual 用summary。函数可以得

17、到回归模型的详细结果，包括系数和上下限。 # 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.560180.48702 0.60199coef()函数得到的系数为向量形式，第一个元素为常数项的系数，第二个及以后为自变量的系数。summary函数se参数的说明：A. se ="rank

18、" ：按照Koenker(1994)的排秩方法计算得到的置信区间，默认残差为独立同分布。注意的是，上下限是不对称的。Coefficients: coefficients lower bd upper bd (Intercept) 81.48225 53.25915 114.01156 income 0.560180.48702 0.60199B. se=" iid " ：假设残差为独立同分布，用KB ( 1978)的方法计算得到近似的协方差矩阵。Coefficients:Value Std. Error t value Pr(>|t|)(Intercept)

19、 81.48225 13.23908 6.15468 0.00000 income 0.56018 0.01192 46.99766 0.00000C. se ="nid” ：表示才照Huber方法逼近得到的估计量。Coefficients:Value Std. Error t value Pr(>|t|)(Intercept) 81.48225 19.25066 4.23270 0.00003income 0.56018 0.02828 19.81032 0.00000D. se=“ ker” :采用 Powell(1990)的核估计方法。Coefficients:Value

20、 Std. Error t value Pr(>|t|)(Intercept) 81.48225 30.21532 2.69672 0.00751 income 0.56018 0.03732 15.01139 0.00000E. 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

21、.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)结果:(krtmefit)图2.1 99个分位点的系数估计值(krtmeiitjOutUI图2.2 5个分

22、位点的系数估计值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),

23、 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)1000200030W 4HH 5DWt iuuxl

24、nld InoMne图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$sol1,#每个分位点的tau值qs.poor =

25、 c( c(1,x.poor) %*% z$sol4:5, ) # 10%分位点的收入的消费估计值qs.rich = c( c(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(q

26、s.poor1,qs.poor), do.points=F,add=T)plot(stepfun(ps,c(qs.poor1,qs.rich), do.points=F,add=T, 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,

27、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")吕Q

28、图2.4 10%分位点和90粉位点之间的比较3-上图表示收入（income）为10%H立点处（poor,穷人）和90%H立点处（rich ,富人）的食品支出的比较。从左图可以发现，对于穷人而言，在不同分位 点估计的食品消费差别不大。而对于富人而言，在不同分位点对食品消费的差别比较大。右图反应了穷人和富人的食品消费分布曲线。穷人的食品消费集中于400左右，比较陡峭；而富人的消费支出集中于 800到1200之间，比较分散。（四）模型比较#比较不同分位点下，收入对食品支出的影响机制是否相同fit1 = rq(foodexp income, tau = 0.25) fit2 = rq(foodexp

29、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.75 Df Resid Df F value Pr(>F)1 2703 15.557 2.449e-07 *Signif. codes: 0 '* ' 0.001 '*&

30、#39; 0.01 '*' 0.05 '.'0.1 ' ' 1其中P值远小于0.05 ,故不同分位点下收入对食品支出的影响机制不同。(五)残差形态的检验也可以理解为是比较不同分位点的模型之间的关系。主要有两种模型形式：(1)位置漂移模型：不同分位点的估计结果之间的 斜率相同或近似，只是截距不同；表现为不同分位点下的拟合曲线是平行的。(2)位置-尺度漂移模型：不同分位点的估计结果之间的 斜率和截距都不同；表现为不同分位点下的拟合曲线不是平行的。# 残差形态的检验source("C:/Program Files/R/R-2.15.0/lib

31、rary/quantreg/doc/gasprice.R")x = gaspricen = length(x)p = 5X = cbind(x(p-1):(n-1),x(p-2):(n-2),x(p-3):(n-3),x(p-4):(n-4)y = xp:n# 位置漂移模型的检验T1 = KhmaladzeTest(yX, taus = -1, nullH="location")T2 = KhmaladzeTest(yX, taus = 10:290/300,nullH="location", se="ker")# 位置尺度

32、漂移模型的检验T3 = KhmaladzeTest(yX, taus = -1, nullH="location-scale")T4 = KhmaladzeTest(yX, taus = 10:290/300,nullH="location-scale", se="ker")结果：运行T1,可以查看其检验结果。其中nullH表示原假设为"location ”,即原假设为位置漂移模型。Tn表示模型整体的检验，统计量为 4.8。THn是对 每个自变量的检验。比较T1和T3的结果(T3的原假设为"位置尺度漂移模型”)，T

33、1的统计量大于T3的统计量，可见相对而言，拒绝"位置漂移模型”的 概率更大，故相对而言“位置尺度漂移模型”更加合适一些。> T1$nullH1 "location" $Tn1 4.803762 $THnX1 X2 X3 X41.0003199 0.5321693 0.5020834 0.8926828attr(,"class")1 "KhmaladzeTest"> T3 $nullH 1 "location-scale" $Tn1 2.705583 $THnX1 X2 X3 X41.21028

34、99 0.6931785 0.5045163 0.8957127attr(,"class")1 "KhmaladzeTest"(六)非线性分位数回归这里的非线性函数为Frank copula函数。# Demo of nonlinear quantile regression model based on Frank copulavFrank <- function(x, df, delta, u) #某个非线性过程，得到的是0,1的值-log(1-(1-exp(-delta)/(1+exp(-delta*pt(x,df)*(1/u)-1)/delt

35、a#非线性模型FrankModel <- function(x, delta, mu,sigma, df, tau) z <- qt(vFrank(x, df, delta, u = tau), df) mu + sigma*zn <- 200 # 样本量df <- 8 # 自由度delta <- 8 #初始参数生成基于T分布的随机数set.seed(1989) x <- sort(rt(n,df) #基于x生成理论上的非参数对应对应的T分布统计量散点图基本数据集v <- vFrank(x, df, delta, u = runif(n) #值y &l

36、t;- qt(v, df)# vwindows(5,5)plot(x, y, pch="o", col="blue", cex = .25) #Dat <- data.frame(x = x, y = y) #us <- c(.25,.5,.75)for(i in 1:length(us)v <- vFrank(x, df, delta, u = usi)lines(x, qt(v,df) # v 为概率，计算每个概率对应的T分布统计量cfMat <- matrix(0, 3, length(us)+1) #初始矩阵，用于保存结果

37、的系数for(i in 1:length(us) tau <- usicat("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

38、, newdata=x), lty=2, col="red") #绘制预测曲线cfMati,1 <- tau #保存分位点的值cfMati,2:4 <- coef(fit) #保存系数到cfMat矩阵的第i行cat("n")#如果前面把每步的结果显示出来，则每次的结果之间添加换行符colnames(cfMat) <- c("分位点",names(coef(fit) #给保存系数的矩阵添加列名cfMat结果:图2.5非参数散点图、真实T分布和拟合结果的比较拟合结果：(过程略)> cfMat分位点 delta mu

39、 sigma1, 0.25 14.87165 -0.20530041 0.91346572, 0.50 16.25362 0.03232525 0.96382093, 0.75 12.09836 0.11998614 0.9423476(七)半参数和非参数分位数回归非参数分位数回归在局部多项式的框架下操作起来更加方便。可以基于以下函数。# 2-7-1半参数模型-fit<-rq(ybs(x,df=5)+z,tau=.33)# 其中bs()表示按b-spline 的非参数拟合# 2-7-2非参数方法lprq <-function(x,y,h,m=50,tau=0.5)#这是自定义的一个

40、非参数计算函数，在其他数据下同样可以使用xx<-seq(min(x),max(x),length=m) # m 个监测点fv<-xxdv<-xx for(i in 1:length(xx)z<-x-xxiwx<-dnorm(z/h) #核函数为正态分布，dnorm计算标准正态分布的密度值r<-rq(yz,weights=wx,tau=tau, #上面计算得到的密度值为权重ci=FALSE)fvi<-r$coef1dvi<-r$coef2list(xx=xx,fv=fv,dv=dv) # 输出结果library(MASS)data(mcycle)a

41、ttach(mcycle)windows(5,5) #非参数的结果一般是通过画图查看的plot(times,accel,xlab="milliseconds",ylab="acceleration")hs<-c(1,2,3,4) #选择不同窗宽进行估计for(i in hs)h=hsifit<-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"

42、,"h=4"),lty=1:length(hs)结果：10 M 3040 SO图2.6基于正态分布为核函数的非参数拟合结果# 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(xhoppers,yho

43、ppers,pch="h",col="red")points(xspecials,yspecials,pch="s",col="blue")others<-(!hoppers&!specials)points(xothers,yothers,col="black",cex=0.75)fk-rqss(yqss(x,lambda=1),tau=0.9) #关键的拟合步骤plotfit,add=T) # add=T表示在上图中添加这里的曲线图2.7 90%分位点下的附加分位数回归拟合结果

44、(A)分位数回归的分解# MM200盼位数分解的函数，改变参数可直接使用MM2005 = function(formu,taus, data, group, pic=F)# furmu 为方程，如 foodexpincome# taus 为不同的分位数# data 总的数据集# group分组指标，是一个向量，用于按行区分data# pic 是否画图，如果分位数比较多，建议不画图engel1 = datagroup=1,engel2 = datagroup=2,# 开始进行分解fita = summary( rq(formu, tau = taus, data = engell )fitb =

45、 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)!=formu2 ) %*%fitai$coef,1yb = cbind(1,

46、engel2,names(engel2)!=formu2 ) %*%fitbi$coef,1# 这里以group=1为基准模型，用group=2的数据计算反常规模型拟合值ystar = cbind(1,engel2,names(engel2)!=formu2 ) %*%fitai$coef,1ya = mean(ya)yb = mean(yb)ystar = mean(ystar)tabi,1 = fitai$tautabi,2 = yb - yatabi,3 = yb - ystar #回报影响，数据相同，模型不同：模型机制的不同所产生的差异tabi,4 = ystar - ya #变量影响

47、，数据不同，模型相同：样本点不同产生的差异#画图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( fitai, col=2 )abline( fitbi, col=3 )detach(engel)#输出结果tab#下面是用一个样本数据做的测试data

48、(engel)group = c(rep(1,100),rep(2,135) # 取前 100 个为第一组，后 135 个第二组taus = c(0.05,0.25,0.5,0.75,0.95) # 需要考察的不同分位点MM2005(foodexpincome, taus, data = engel, group=group, pic=F) #参数说明，见说明： 自编函数MM2005勺参数说明：函数调用形式 MM2005 (formu, taus, data, group, pic=F),其中# furmu 为回归方程，如foodexpincome;# taus 为不同的分位数，如 taus=

49、c(0.05,0.5,0.95);# data 总的数据集，如上例中的engel数据集；# group分组指标，是一个向量，用于按行区分 data,第一组为1,第二组为2;目前仅能分两组；# pic逻辑参数：是否画图。如果分位数比较多，建议不画图；默认不画图，pic=F;如果想画图，则添加参数pic=T o最终结果：> MM2005(foodexpincome, taus, data = engel, group=group, pic=F) 分位数 总差异回报影响变量影响0.05 -30.452061 -72.35939 41.907330.25 -2.017317 -46.20125

50、44.183940.50 30.941212 -23.24042 54.181630.75 43.729025 -15.76283 59.491850.95 52.778985 -11.29932 64.07830第三节一个案例这里使用某年份流动人口中个人消费和个人收入的数据作为案例，进行整体性分析。(一)初始化程序基础数据是dta格式(stata软件的数据格式)，所以要用foreign包里的read.dta函数读取。这部分的主要内容包括：1、修改工作目录至程序和数据所在文件夹，注意不能使用""，必须使用"/"或""；2、清除其他变量

51、、清理内存；3、加载分位数回归所需的包quantreg和导入数据所需的foreign包；4、加载4-functions.R 文件中的R程序，这里有三个有用的函数；5、读取案例数据20121218.dta ;dat26、初始化处理：去掉含有缺失值的行，检查样本个数，检查列名，整理列名，进行对数变换并生成新的数据集# 一个案例setwd("F:/实践2/2012_分位数回归等两个任务/分位数回归方法夹")#设置工作目录setwd(" " )rm(list=ls()gc()library(quantreg)library(foreign)source(&quo

52、t;4-functions.R") # 将 4-functions.R 中的函数放入内存中dat = read.dta("20121218.dta")dat = na.omit(dat) #去掉有缺失值的行dim(dat)# 这个数据包含10000个样本，3个指标names(dat)# 1 "urtype” "q413" "q410a"# q413每月家庭总收入# q410a家庭食品支出# urtype 样本类型，1农村或2城镇# 这里考察上个月的收入与食品支出的关系，并分性别进行讨论# 把变量的名称改一下name

53、s(dat) = c("urtype","income","foodexp")# 进行log变换dat2 = datdat$income>0 & dat$foodexp>0,dat2,"lnincome" = log( dat2,"income")dat2,"lnfoodexp" = log( dat2,"foodexp")dat2 = na.omit(dat2) # 去掉log变换以后的缺失值行dim(dat2)(二)画散点图画散点图

54、的主要目的是检查数据的分布是否需要进行对数变换，另外查看一下使用什么样的模型更合适。本案例中进行对数变换是要合适一些。# 散点图，观察是否需要log变换attach(dat)windows(5,5)# 不变换plot(income,foodexp,cex=0.5,main=" 散点图")detach(dat)# 变换attach(dat2)windows(5,5)plot(lnincome,lnfoodexp,cexp=0.5,main=" 散点图(对数变换以后)")detach(dat2)结果:敢苴附iiKone数点图(对敷变摸以息J34567 B J

55、mrcorneI图3.1进行对数变换或不对数变换下的散点图(三)分位数回归的参数估计结果调用rq函数进行分位数回归估计。1、这里选取的分位点包括0.05、0.25、0.5、0.75和0.9。分别对整体数据、城镇和农村分别建模。如果数据量较大，建议method参数为"fn"或"pfn", 运算效率会高些。2、用summary。函数可以得到更加详细的统计结果。3、用笔者编写的tabs()函数和tabcoef()函数可以把summary。处理过的结果进行整合。注意,笔者的这两个函数在4-funcitons.R 中，只有加载了这个文件(source( "

56、;4-functions.R ")才能使用这两个函数，#建立基础的分位数回归模型taus = c(0.05,0.25,0.5,0.75,0.95) fitl = rq( Infoodexp Inincome, data=dat2, tau=taus, method="fn")fit2 = rq( lnfoodexp lnincome, data=dat2dat2$urtype=1,tau=taus, method="fn")fit3 = rq( lnfoodexp lnincome, data=dat2dat2$urtype=2,tau=taus, method="fn")#所有数据放在一起是参数结果si = summary(fitl)#男性s2 = summary(fit2)#女性s3 = summary(fit3)#参数结果的直接比较tabs(s1) # 这里的tabs函数是笔