用R語言進(jìn)行分位數(shù)回歸

本文格式為Word版，下載可任意編輯——用R語言進(jìn)行分位數(shù)回歸|t|)

(Intercept)81.4822513.239086.154680.00000income0.560180.0119246.997660.00000C.se=“nid〞:表示依照Huber方法迫近得到的估計(jì)量。Coefficients:

ValueStd.ErrortvaluePr(>|t|)

(Intercept)81.4822519.250664.232700.00003income0.560180.0282819.810320.00000D.se=〞ker〞:采用Powell(1990)的核估計(jì)方法。Coefficients:

ValueStd.ErrortvaluePr(>|t|)

(Intercept)81.4822530.215322.696720.00751income0.560180.0373215.011390.00000E.se=〞boot〞:采用bootstrap方法自助抽樣的方法估計(jì)系數(shù)的誤差標(biāo)準(zhǔn)差。Coefficients:ValueStd.ErrortvaluePr(>|t|)(Intercept)81.4822525.236473.228750.00142income0.560180.0319417.537520.00000（三）不同分位點(diǎn)下的回歸結(jié)果比較1、不同分為點(diǎn)系數(shù)估計(jì)值的比較#不同分位點(diǎn)下的系數(shù)估計(jì)值的比較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)#新建一個(gè)圖形窗口，可以去掉這句plot(fit1)windows(5,5)#新建一個(gè)圖形窗口，可以去掉這句plot(fit2)結(jié)果：圖2.199個(gè)分位點(diǎn)的系數(shù)估計(jì)值圖2.25個(gè)分位點(diǎn)的系數(shù)估計(jì)值2、不同分位點(diǎn)擬合曲線的比較#散點(diǎn)圖attach(engel)#開啟engel數(shù)據(jù)集，直接運(yùn)行其中的列名，就可以調(diào)用相應(yīng)列plot(income,foodexp,cex=0.25,type=\#畫圖，說明①xlab=\points(income,foodexp,cex=0.5,col=\#添加點(diǎn)，點(diǎn)的大小為0.5abline(rq(foodexp~income,tau=0.5),col=\#畫中位數(shù)回歸的擬合直線，顏色藍(lán)abline(lm(foodexp~income),lty=2,col=\#畫普通最小二乘法擬合直線，顏色紅taus=c(0.05,0.1,0.25,0.75,0.9,0.95)for(iin1:length(taus)){#繪制不同分位點(diǎn)下的擬合直線，顏色為灰色abline(rq(foodexp~income,tau=taus[i]),col=\}detach(engel)圖2.3不同分位點(diǎn)下的分位數(shù)回歸擬合結(jié)果比較3、窮人和富人的消費(fèi)分布比較#比較窮人（收入在10%分位點(diǎn)的那個(gè)人）和富人（收入在90%分位點(diǎn)的那個(gè)人）的估計(jì)結(jié)果#rq函數(shù)中，tau不在[0,1]時(shí)，表示按最細(xì)的分位點(diǎn)劃分方式得到分位點(diǎn)序列z=rq(foodexp~income,tau=-1)z$sol#這里包含了每個(gè)分位點(diǎn)下的系數(shù)估計(jì)結(jié)果x.poor=quantile(income,0.1)#10%分位點(diǎn)的收入x.rich=quantile(income,0.9)#90%分位點(diǎn)的收入ps=z$sol[1,]#每個(gè)分位點(diǎn)的tau值qs.poor=c(c(1,x.poor)%*%z$sol[4:5,])#10%分位點(diǎn)的收入的消費(fèi)估計(jì)值qs.rich=c(c(1,x.rich)%*%z$sol[4:5,])#90%分位點(diǎn)的收入的消費(fèi)估計(jì)值windows(10,5)par(mfrow=c(1,2))#把繪圖區(qū)域劃分為一行兩列plot(c(ps,ps),c(qs.poor,qs.rich),type=\#type=〞n〞表示初始化圖形區(qū)域，但不畫圖xlab=expression(tau),ylab=\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=\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=\lines(qs.rich,ar$dens,col=\lines(qs.poor,ap$dens,col=\legend(\col=c(\圖2.410%分位點(diǎn)和90%分位點(diǎn)之間的比較上圖表示收入（income）為10%分位點(diǎn)處（poor，窮人）和90%分位點(diǎn)處（rich，富人）的食品支出的比較。從左圖可以發(fā)現(xiàn)，對(duì)于窮人而言，在不同分位點(diǎn)估計(jì)的食品消費(fèi)區(qū)別不大。而對(duì)于富人而言，在不同分位點(diǎn)對(duì)食品消費(fèi)的區(qū)別比較大。右圖反應(yīng)了窮人和富人的食品消費(fèi)分布曲線。窮人的食品消費(fèi)集中于400左右，比較陡峭；而富人的消費(fèi)支出集中于800到1200之間，比較分散。（四）模型比較#比較不同分位點(diǎn)下，收入對(duì)食品支出的影響機(jī)制是否一致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)結(jié)果：QuantileRegressionAnalysisofDevianceTableModel:foodexp~incomeJointTestofEqualityofSlopes:tauin{0.250.50.75}DfResidDfFvaluePr(>F)1270315.5572.449e-07***Signif.codes:0?***?0.001?**?0.01?*?0.05?.?0.1??1其中P值遠(yuǎn)小于0.05，故不同分位點(diǎn)下收入對(duì)食品支出的影響機(jī)制不同。（五）殘差形態(tài)的檢驗(yàn)也可以理解為是比較不同分位點(diǎn)的模型之間的關(guān)系。主要有兩種模型形式：

（1）位置漂移模型：不同分位點(diǎn)的估計(jì)結(jié)果之間的斜率一致或近似，只是截距不同；表現(xiàn)為不同分位點(diǎn)下的擬合曲線是平行的。

（2）位置-尺度漂移模型：不同分位點(diǎn)的估計(jì)結(jié)果之間的斜率和截距都不同；表現(xiàn)為不同分位點(diǎn)下的擬合曲線不是平行的。

#殘差形態(tài)的檢驗(yàn)

source(\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=x[p:n]

#位置漂移模型的檢驗(yàn)

T1=KhmaladzeTest(y~X,taus=-1,nullH=\T2=KhmaladzeTest(y~X,taus=10:290/300,nullH=\#位置尺度漂移模型的檢驗(yàn)

T3=KhmaladzeTest(y~X,taus=-1,nullH=\T4=KhmaladzeTest(y~X,taus=10:290/300,nullH=\se=\

結(jié)果：運(yùn)行T1，可以查看其檢驗(yàn)結(jié)果。其中nullH表示原假設(shè)為“l(fā)ocation〞，即原假設(shè)為位置漂移模型。Tn表示模型整體的檢驗(yàn)，統(tǒng)計(jì)量為4.8。THn是對(duì)每個(gè)

自變量的檢驗(yàn)。比較T1和T3的結(jié)果（T3的原假設(shè)為“位置尺度漂移模型〞），T1的統(tǒng)計(jì)量大于T3的統(tǒng)計(jì)量，可見相對(duì)而言，拒絕“位置漂移模型〞的概率更大，

故相對(duì)而言“位置尺度漂移模型〞更加適合一些。

>T1$nullH[1]\$Tn[1]4.803762$THn

X1X2X3X4

1.00031990.53216930.50208340.8926828attr(,\[1]\>T3$nullH

[1]\$Tn[1]2.705583$THn

X1X2X3X4

1.21028990.69317850.50451630.8957127

attr(,\[1]\

（六）非線性分位數(shù)回歸這里的非線性函數(shù)為Frankcopula函數(shù)。##DemoofnonlinearquantileregressionmodelbasedonFrankcopulavFrankcfMat分位點(diǎn)deltamusigma[1,]0.2514.87165-0.205300410.9134657[2,]0.5016.253620.032325250.9638209[3,]0.7512.098360.119986140.9423476（七）半?yún)?shù)和非參數(shù)分位數(shù)回歸非參數(shù)分位數(shù)回歸在局部多項(xiàng)式的框架下操作起來更加便利?？梢曰谝韵潞瘮?shù)。#2-7-1半?yún)?shù)模型fit

nullH=\T2=KhmaladzeTest(lnfoodexp~lnincome,taus=seq(.1,.9,by=.1),nullH=\T1T2detach(dat2)結(jié)果：根據(jù)以下結(jié)果，以位置漂移模型為原假設(shè)的統(tǒng)計(jì)量為2.84，以位置漂移模型為原假設(shè)的統(tǒng)計(jì)量為3.29,。相對(duì)而言，更簡(jiǎn)單拒絕位置漂移模型，故依照位置漂移模型相對(duì)適合一些，即不同分位點(diǎn)的斜率基本上可以認(rèn)為是一致的，而截距不同。>T1$nullH[1]\$Tn[1]2.840216$THn[1]2.840216attr(,\[1]\>T2$nullH[1]\$Tn[1]3.2899$THn[1]3.2899attr(,\[1]\（五）分位數(shù)回歸的圖形結(jié)果#畫圖#參數(shù)比較windows(5,5)plot(s1)windows(5,5)plot(s2)windows(5,5)plot(s3)#散點(diǎn)圖#-總體windows(5,5)attach(dat2)plot(lnincome,lnfoodexp,cex=0.25,type=\xlab=\main=\points(lnincome,lnfoodexp,pch=\abline(rq(lnfoodexp~lnincome,tau=0.5,method=\abline(lm(lnfoodexp~lnincome),lty=2,col=\taus=c(0.05,0.1,0.25,0.75,0.9,0.95)for(iin1:length(taus)){abline(rq(lnfoodexp~lnincome,tau=taus[i],method=\col=\}legend(\\detach(dat2)#-城鎮(zhèn)windows(5,5)dat2a=dat2[dat2$urtype==1,]attach(dat2a)plot(lnincome,lnfoodexp,cex=0.25,type=\xlab=\main=\points(lnincome,lnfoodexp,pch=\abline(rq(lnfoodexp~lnincome,tau=0.5,method=\abline(lm(lnfoodexp~lnincome),lty=2,col=\taus=c(0.05,0.1,0.25,0.75,0.9,0.95)for(iin1:length(taus)){abline(rq(lnfoodexp~lnincome,tau=taus[i],method=\col=\}legend(\\detach(dat2a)#-農(nóng)村windows(5,5)dat2b=dat2[dat2$urtype==2,]attach(dat2b)plot(lnincome,lnfoodexp,cex=0.25,type=\xlab=\main=\points(lnincome,lnfoodexp,pch=\abline(rq(lnfoodexp~lnincome,tau=0.5,method=\abline(lm(lnfoodexp~lnincome),lty=2,col=\taus=c(0.05,0.1,0.25,0.75,0.9,0.95)for(iin1:length(taus)){abline(rq(lnfoodexp~lnincome,tau=taus[i],method=\col=\}legend(\\detach(dat2b)結(jié)果：圖3.2.a整體數(shù)據(jù)的分位數(shù)回歸的系數(shù)走勢(shì)圖3.2.b城鎮(zhèn)數(shù)據(jù)的分位數(shù)回歸的系數(shù)走勢(shì)

圖3.2.c農(nóng)村數(shù)據(jù)的分位數(shù)回歸的系數(shù)走勢(shì)

圖3.3.a全部數(shù)據(jù)的分位數(shù)回歸的散點(diǎn)圖

圖3.3.b城鎮(zhèn)數(shù)據(jù)的分位數(shù)回歸的散點(diǎn)圖圖3.3.c農(nóng)村數(shù)據(jù)的分位數(shù)回歸的散點(diǎn)圖（六）分位數(shù)分解依照MM2023方法進(jìn)行分解。注意的是：1、group參數(shù)必需是包含1和2的向量，其中反事實(shí)模型的擬合值是采用1類的模型參數(shù)和2類的數(shù)據(jù)計(jì)算得到。相應(yīng)地，在本案例中作為參照的擬合結(jié)果是假定農(nóng)村滾動(dòng)人口依照城鎮(zhèn)的決定模型，會(huì)得到什么樣的食品支出。2、data參數(shù)中只能包含模型中出現(xiàn)的變量。譬如這里使用的變量只包括lnfoodexp和lnincome，所以只能保存這兩個(gè)變量在data中。#依照城鄉(xiāng)進(jìn)行分解group=as.character(dat2$urtype)#必需保證group是包含1和2的向量group[group==\group[group==\taus=c(0.05,0.25,0.5,0.75,0.95)data=dat2[,c(\res=MM2023(lnfoodexp~lnincome,taus,data,group,pic=F,m.rq=\#分位數(shù)總差異回報(bào)影響變量影響#0.05-0.0510898590.06520406-0.11629392#0.25-0.0565060900.04647729-0.10298338#0.50-0.0146327740.08108505-0.09571782#0.750.0096296910.10826533-0.09863564#0.950.0287912460.10291667-0.07412543write.csv(res,\#將結(jié)果保存為csv表格，便利整理結(jié)果：表3-2分位數(shù)分解：按城鄉(xiāng)的分解分位數(shù)0.050.250.500.750.95

總差異-0.22-0.18-0.23-0.27-0.41

回報(bào)影響-0.04-0.01-0.06-0.10-0.23

變量影響-0.18-0.17-0.17-0.17-0.17

以上結(jié)果顯示，無論對(duì)于高收入家庭還是低收入家庭，個(gè)體對(duì)食品支出的影響都差不多。但對(duì)于高收入家庭