用R语言进行分位数回归.docx
《用R语言进行分位数回归.docx》由会员分享，可在线阅读，更多相关《用R语言进行分位数回归.docx（26页珍藏版）》请在冰豆网上搜索。
用R语言进行分位数回归
用R语言进行分位数回归：
基础篇
第一节 分位数回归介绍
[2][3]
[4]
[5][6]
[7]
[8]
[9]
[10]
[11]
第二节 用R语言进行分位数回归
（“quantreg”） #保持联网的情况下安装包
library（“quantreg”） #加载包
（） #进入R帮助首页
help（rq） #获取rq函数的帮助，也可以写成：
rq
example（rq） #显示分位数回归函数rq（）的一个简单示例代码
data（engel） #加载quantreg包自带的数据集，见说明①
fit1=rq（foodexp~income,tau=,data=engel） #进行分位数回归，见说明②
fit1 #直接显示分位数回归的模型和系数，见说明③
summary（fit1） #得到更加详细的显示结果，见说明④
r1=resid（fit1） #得到残差序列，并赋值为变量r1
c1=coef（fit1） #得到模型的系数，并赋值给变量c1，见说明⑤
summary（fit1,se=“nid”） #通过设置参数se，可以得到系数的假设检验，说明⑥
Frisch–Newton内点算法；C.“pfn”，针对特别大数据，使用经过预处理的Frisch–Newton逼近方法；D.“fnc”，针对被拟合系数特殊的线性不等式约束情况；E.“lasso”和“scad”，基于特定惩罚函数的平滑算法进行拟合。
#不同分位点下的系数估计值的比较
fit1=summary（rq（foodexp~income,tau=2:
98/100））
fit2=summary（rq（foodexp~income,tau=c,,,,））
windows（5,5） #新建一个图形窗口，可以去掉这句
plot（fit1）
windows（5,5） #新建一个图形窗口，可以去掉这句
plot（fit2）
#散点图
attach（engel） #打开engel数据集，直接运行其中的列名，就可以调用相应列
plot（income,foodexp,cex=,type="n", #画图，说明①
xlab="HouseholdIncome",ylab="FoodExpenditure"）
points（income,foodexp,cex=,col="blue"） #添加点，点的大小为
abline（rq（foodexp~income,tau=,col="blue"） #画中位数回归的拟合直线，颜色蓝
abline（lm（foodexp~income）,lty=2,col="red"） #画普通最小二乘法拟合直线，颜色红
taus=c,,,,,
for（iin1:
length（taus））{ #绘制不同分位点下的拟合直线，颜色为灰色
abline（rq（foodexp~income,tau=taus[i]）,col="gray"）
}
detach（engel）
#比较穷人（收入在10%分位点的那个人）和富人（收入在90%分位点的那个人）的估计结果
#rq函数中，tau不在[0,1]时，表示按最细的分位点划分方式得到分位点序列
z=rq（foodexp~income,tau=-1）
z$sol #这里包含了每个分位点下的系数估计结果
=quantile（income,#10%分位点的收入
=quantile（income,#90%分位点的收入
ps=z$sol[1,] #每个分位点的tau值
=c（c（1,%*%z$sol[4:
5,]） #10%分位点的收入的消费估计值
=c（c（1,%*%z$sol[4:
5,]） #90%分位点的收入的消费估计值
windows（10,5）
par（mfrow=c（1,2）） #把绘图区域划分为一行两列
plot（c（ps,ps）,c,,type="n", #type=”n”表示初始化图形区域，但不画图
xlab=expression（tau）,ylab="quantile"）
plot（stepfun（ps,c[1],）,=F,
add=T）
plot（stepfun（ps,c[1],）,=F,
add=T,="gray",="gray"）
=（c（0,diff（ps））+c（diff（ps）,0））/2
ap=akj,z=,p=
ar=akj,z=,p=
plot（c,,c（ap$dens,ar$dens）,
type="n",xlab="FoodExpenditure",ylab="Density"）
lines,ar$dens,col="gray"）
lines,ap$dens,col="black"）
legend（"topright",c（"poor","rich"）,lty=c（1,1）,
col=c（"black","gray"））
#比较不同分位点下，收入对食品支出的影响机制是否相同
fit1=rq（foodexp~income,tau=
fit2=rq（foodexp~income,tau=
fit3=rq（foodexp~income,tau=
anova（fit1,fit2,fit3）
#残差形态的检验
source（"C:
/ProgramFiles/R/"）
x=gasprice
n=length（x）
p=5
X=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]
#位置漂移模型的检验
T1=KhmaladzeTest（y~X,taus=-1,nullH="location"）
T2=KhmaladzeTest（y~X,taus=10:
290/300,
nullH="location",se="ker"）
#位置尺度漂移模型的检验
T3=KhmaladzeTest（y~X,taus=-1,nullH="location-scale"）
T4=KhmaladzeTest（y~X,taus=10:
290/300,
nullH="location-scale",se="ker"）
##DemoofnonlinearquantileregressionmodelbasedonFrankcopula
vFrank<-function（x,df,delta,u） #某个非线性过程，得到的是[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*z
}
n<-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,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=us[i]）
lines（x,qt（v,df）） #v为概率，计算每个概率对应的T分布统计量
}
cfMat<-matrix（0,3,length（us）+1）#初始矩阵，用于保存结果的系数
for（iin1:
length（us））{
tau<-us[i]
cat（"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,newdata=x）,lty=2,col="red"） #绘制预测曲线
cfMat[i,1]<-tau #保存分位点的值
cfMat[i,2:
4]<-coef（fit） #保存系数到cfMat矩阵的第i行
cat（"\n"） #如果前面把每步的结果显示出来，则每次的结果之间添加换行符
}
colnames（cfMat）<-c（"分位点",names（coef（fit））） #给保存系数的矩阵添加列名
cfMat
#2-7-1半参数模型----
fit<-rq（y~bs（x,df=5）+z,tau=.33）
#其中bs（）表示按b-spline的非参数拟合
#2-7-2非参数方法
lprq<-function（x,y,h,m=50,tau={
#这是自定义的一个非参数计算函数，在其他数据下同样可以使用
xx<-seq（min（x）,max（x）,length=m） #m个监测点
fv<-xx
dv<-xx
for（iin1:
length（xx））{
z<-x-xx[i]
wx<-dnorm（z/h） #核函数为正态分布,dnorm计算标准正态分布的密度值
r<-rq（y~z,weights=wx,tau=tau, #上面计算得到的密度值为权重
ci=FALSE）
fv[i]<-r$coef[1]
dv[i]<-r$coef[2]
}
list（xx=xx,fv=fv,dv=dv） #输出结果
}
library（MASS）
data（mcycle）
attach（mcycle）
windows（5,5） #非参数的结果一般是通过画图查看的
plot（times,accel,xlab="milliseconds",ylab="acceleration"）
hs<-c（1,2,3,4） #选择不同窗宽进行估计
for（iinhs）{
h=hs[i]
fit<-lprq（times,accel,h=h,tau= #关键拟合函数
lines（fit$xx,fit$fv,lty=i）
}
legend（45,-70,c（"h=1","h=2","h=3","h=4"）,
lty=1:
length（hs））
#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（x[hoppers],y[hoppers],pch="h",col="red"）
points（x[specials],y[specials],pch="s",col="blue"）
others<-（!
hoppers&!
specials）
points（x[others],y[others],col="black",cex=
fit<-rqss（y~qss（x,lambda=1）,tau= #关键的拟合步骤
plot（fit,add=T） #add=T表示在上图中添加这里的曲线
#MM2005分位数分解的函数，改变参数可直接使用
MM2005=function（formu,taus,data,group,pic=F）{
#furmu为方程，如foodexp~income
#taus为不同的分位数
#data总的数据集
#group分组指标，是一个向量，用于按行区分data
#pic是否画图，如果分位数比较多，建议不画图
engel1=data[group==1,]
engel2=data[group==2,]
#开始进行分解
fita=summary（rq（formu,tau=taus,data=engel1））
fitb=summary（rq（formu,tau=taus,data=engel2））
tab=matrix（0,length（taus）,4）
colnames（tab）=c（"分位数","总差异","回报影响","变量影响"）
rownames（tab）=rep（"",dim（tab）[1]）
for（iin1:
length（taus））{
ya=cbind（1,engel1[,names（engel1）!
=formu[[2]]]）%*%fita[[i]]$coef[,1]
yb=cbind（1,engel2[,names（engel2）!
=formu[[2]]]）%*%fitb[[i]]$coef[,1]
#这里以group==1为基准模型，用group==2的数据计算反常规模型拟合值
ystar=cbind（1,engel2[,names（engel2）!
=formu[[2]]]）%*%fita[[i]]$coef[,1]
ya=mean（ya）
yb=mean（yb）
ystar=mean（ystar）
tab[i,1]=fita[[i]]$tau
tab[i,2]=yb-ya
tab[i,3]=yb-ystar#回报影响，数据相同，模型不同：
模型机制的不同所产生的差异
tab[i,4]=ystar-ya#变量影响，数据不同，模型相同：
样本点不同产生的差异
}
#画图
if（pic）{
attach（engel）
windows（5,5）
plot（income,foodexp,cex=,type="n",main="两组分位数回归结果比较"）
points（engel1,cex=,col=2）
points（engel2,cex=,col=3）
for（iin1:
length（taus））{
abline（fita[[i]],col=2）
abline（fitb[[i]],col=3）
}
detach（engel）
}
#输出结果
tab
}
#下面是用一个样本数据做的测试
data（engel）
group=c（rep（1,100）,rep（2,135）） #取前100个为第一组，后135个第二组
taus=c,,,, #需要考察的不同分位点
MM2005（foodexp~income,taus,data=engel,group=group,pic=F）#参数说明，见①
第三节 一个案例
#一个案例
setwd（"F:
/实践2/2012_分位数回归等两个任务/分位数回归方法夹"）#设置工作目录
setwd（””）
rm（list=ls（））
gc（）
library（quantreg）
library（foreign）
source（""）#将中的函数放入内存中
dat=（".dta"）
dat=（dat） #去掉有缺失值的行
dim（dat）
#这个数据包含10000个样本，3个指标
names（dat）
#[1]"urtype""q413" "q410a"
#q413每月家庭总收入
#q410a家庭食品支出
#urtype样本类型，1农村或2城镇
#这里考察上个月的收入与食品支出的关系，并分性别进行讨论
#把变量的名称改一下
names（dat）=c（"urtype","income","foodexp"）
#进行log变换
dat2=dat[dat$income>0&dat$foodexp>0,]
dat2[,"lnincome"]=log（dat2[,"income"]）
dat2[,"lnfoodexp"]=log（dat2[,"foodexp"]）
dat2=（dat2） #去掉log变换以后的缺失值行
dim（dat2）
#散点图，观察是否需要log变换
attach（dat）
windows（5,5）
#不变换
plot（income,foodexp,cex=,main="散点图"）
detach（dat）
#变换
attach（dat2）
windows（5,5）
plot（lnincome,lnfoodexp,cexp=,main="散点图（对数变换以后）"）
detach（dat2）
#建立基础的分位数回归模型
taus=c,,,,
fit1=rq（lnfoodexp~lnincome,data=dat2,tau=taus,method="fn"）
fit2=rq（lnfoodexp~lnincome,data=dat2[dat2$urtype==1,],
tau=taus,method="fn"）
fit3=rq（lnfoodexp~lnincome,data=dat2[dat2$urtype==2,],
tau=taus,method="fn"）
#所有数据放在一起是参数结果
s1=summary（fit1）
#男性
s2=summary（fit2）
#女性
s3=summary（fit3）
#参数结果的直接比较
tabs（s1） #这里的tabs函数是笔者自己写的，方便显示结果
tabs（s2）
tabs（s3）
#具体的系数估计值及假设检验
tabcoef（s1） #这里的tabcoef函数是笔者自己写的，方便显示结果
tabcoef（s2）
tabcoef（s3）
Quantiles
（Intercept）
lnincome
Quantiles
（Intercept）
lnincome
Quantiles
（Intercept）
lnincome
Quantiles
Value
Std.Error
tvalue
Pr（>|t|）
Quantiles
Value
Std.Error
tvalue
Pr（>|t|）
Quantiles
Value
Std.Error
tvalue
Pr（>|t|）
#参数结果的检验
attach（dat2）
T1=KhmaladzeTest（lnfoodexp~lnincome,taus=seq（.1,.9,by=.1）,
nullH="location",se="ker"）
T2=KhmaladzeTest（lnfoo
升级会员 