《用R语言进行分位数回归》由会员分享，可在线阅读，更多相关《用R语言进行分位数回归（32页珍藏版）》请在金锄头文库上搜索。
1、用R语言进行分位数回归：基础篇91011第二节用R语言进行分位数回归(“quantreg”)#保持联网的情况下安装包library(“quantreg”)#加载包()#进入R帮助首页help(rq)#获取rq函数的帮助，也可以写成：rqexample(rq)#显示分位数回归函数rq()的一个简单示例代码data(engel)#加载quantreg包自带的数据集，见说明fitl二rq(foodexpincome,tau二，data=engel)#进行分位数回归，见说明fitl#直接显示分位数回归的模型和系数，见说明summary(fitl)#得到更加详细的显示结果，见说明rl=resid(fit
2、l)#得到残差序列，并赋值为变量rlcl=coef(fitl)#得到模型的系数，并赋值给变量cl，见说明summary(fitl,se=“nid”)#通过设置参数se，可以得到系数的假设检验，说明Frisch-Newton内点算法；C.“pfn”,针对特别大数据，使用经过预处理的Frisch-Newton逼近方法；D.“fnc”，针对被拟合系数特殊的线性不等式约束情况；E.“lasso”和“scad”，基于特定惩罚函数的平滑算法进行拟合。#不同分位点下的系数估计值的比较fitl二summary(rq(foodexpincome,tau=2:98/100)fit2=summary(rq(food
3、expincome,tau=c,)windows(5,5)plot(fit1)#新建一个图形窗口，可以去掉这句windows(5,5)plot(fit2)#新建一个图形窗口，可以去掉这句a.a(tataTeptJ0-20-40_Bas#散点图attach(engel)#打开engel数据集，直接运行其中的列名，就可以调用相应列plot(income,foodexp,cex二,type=n,#画图，说明xlab二HouseholdIncome,ylab二FoodExpenditure)points(income,foodexp,cex二，col二blue)#添加点，点的大小为abline(rq(
4、foodexpincome,tau=,col二blue)#画中位数回归的拟合直线，颜色蓝abline(lm(foodexpincome),lty=2,col二red)#画普通最小二乘法拟合直线，颜色红taus=c,for(iin1:length(taus)#绘制不同分位点下的拟合直线，颜色为灰色abline(rq(foodexpincome,tau二tausi),col二gray)detach(engel)10002U003HM40009MHotEeholdhcome#比较穷人（收入在10%分位点的那个人）和富人（收入在90%分位点的那个人）的估计结果#rq函数中，tau不在0,1时，表示按最
5、细的分位点划分方式得到分位点序列z=rq（foodexpincome,tau=-1）z$sol#这里包含了每个分位点下的系数估计结果# 10%分位点的收入# 90%分位点的收入)#每个分位点的tau值#10%分位点的收入的消费估计值#90%分位点的收入的消费估计值#把绘图区域划分为#type=”n”表示初始化图形区域，但ylab=quantile)gray(0,diff(ps)+c(diff(ps),0)/2z=,p=z=,p=c(ap$dens,ar$dens),=quantile(income,=quantile(income,ps=z$sol1,=c(c(1,%*%z$sol4:5,=c
6、(c(1,%*%z$sol4:5,windows(10,5)par(mfrow二c(l,2)一行两列plot(c(ps,ps),c,type二n，不画图xlab二expression(tau),plot(stepfun(ps,c1,),=F,add二T)plot(stepfun(ps,c1,),=F,add=T,二gray,=(cap二akj,ar=akj,plot(c,type=n,xlab二FoodExpenditure,ylab=Density)lines,ar$dens,col二gray)lines,ap$dens,col二black)legend(topright，c(poor，ri
7、ch)，lty二c(l,l),col二c(black，gray)#比较不同分位点下，收入对食品支出的影响机制是否相同fitl二rq(foodexpincome,tau=fit2=rq(foodexpincome,tau=fit3=rq(foodexpincome,tau=anova(fit1,fit2,fit3)#残差形态的检验source(C:/ProgramFiles/R/)x=gaspricen=length(x)p=5X=cbind(x(pl):(nT),x(p2):(n2),x(p3):(n3),x(p4):(n4)y=xp:n#位置漂移模型的检验T1=KhmaladzeTest(y
8、X,taus=1,nullH=location)T2=KhmaladzeTest(yX,taus=10:290/300,nullH=location,se=ker)#位置尺度漂移模型的检验T3=KhmaladzeTest(yX,taus=1,nullH=location-scale)T4=KhmaladzeTest(yX,taus=10:290/300,nullH=location-scale,se=ker)#DemoofnonlinearquantileregressionmodelbasedonFrankcopulavFrank-function(x,df,delta,u)#某个非线性过程
9、，得到的是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),df)mu+sigma*zn-200#样本量df-8#自由度delta-8#初始参数(1989)x-sort(rt(n,df)#生成基于T分布的随机数v-vFrank(x,df,delta,u=runif(n)#基于x生成理论上的非参数对应值y-qt(v,df)#v对应的T分布统计量windows(5
10、,5)plot(x,y,pch=o,col二blue,cex=.25)#散点图Dat-(x=x,y=y)#基本数据集us-c(.25,.5,.75)for(iin1:length(us)v-vFrank(x,df,delta,u=usi)lines(x,qt(v,df)#v为概率，计算每个概率对应的T分布统计量cfMat-matrix(0,3,length(us)+l)#初始矩阵，用于保存结果的系数for(iin1:length(us)tau-usicat(tau=,format(tau),.)fit一nlrq(yFrankModel(x,delta,mu,sigma,df=8,tau=tau
11、),#非参数模型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)#绘制预测曲线cfMati,l-tau#保存分位点的值cfMati,2:4-coef(fit)#保存系数到cfMat矩阵的第i行#2-7-2非参数方法lprq-function(x,y,h,m=50,tau=#这是自定义的一个非参数计算函数，在其他数据下同样可以使用xx-seq(min(x),
12、max(x),length二m)#m个监测点fv-xxdv-xxfor(iin1: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)attach(mcycle)windows(5,5)#非参数的结果一般是通过画图查看的plot(times,accel,xlab=milliseconds
13、,ylab=acceleration)hs-c(l,2,3,4)#选择不同窗宽进行估计for(iinhs)h=hsifit-lprq(times,accel,h二h,tau=#关键拟合函数lines(fit$xx,fit$fv,lty二i)legend(45,-70,c(h=l,h=2,h=3,h=4),lty二l:length(hs)#2-7-3非参数回归的另一个方法#考察最大的跑步速度与体重的关系data(Mammals)attach(Mammals)x-log(weight)#取得自变量的值y-log(speed)#取得因变量的值plot(x,y,xlab二Weightinlog(Kg)
14、,ylab二Speedinlog(Km/hour),type二n)points(xhoppers,yhoppers,pch=h,col=red)points(xspecials,yspecials,pch二s，col二blue)others-(!hoppers&!specials)points(xothers,yothers,col=black,cex=fit-rqss(yqss(x,lambda=1),tau二#关键的拟合步骤plot(fit,add二T)#add=T表示在上图中添加这里的曲线EffectofxIIIIIIIIIIIII11ihll1111111111111111IIIIIIIII1II11111111-20246#MM2005分位数分解的函数，改变参数可直接使用MM2005=function(fo
