分位数回归.docx
《分位数回归.docx》由会员分享，可在线阅读，更多相关《分位数回归.docx（32页珍藏版）》请在冰点文库上搜索。
分位数回归
2、不同分位点拟合曲线的比较
# 散点图
attach（engel） # 打开engel数据集，直接运行其中的列名，就可以调用相应列
plot（income,foodexp,cex=0.25,type="n", # 画图，说明①
xlab="HouseholdIncome",ylab="FoodExpenditure"）
points（income,foodexp,cex=0.5,col="blue"） # 添加点，点的大小为0.5
abline（rq（foodexp~income,tau=0.5）,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（iin1:
length（taus））{ # 绘制不同分位点下的拟合直线，颜色为灰色
abline（rq（foodexp~income,tau=taus[i]）,col="gray"）
}
detach（engel）
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$sol[1,] # 每个分位点的tau值
qs.poor=c（c（1,x.poor）%*%z$sol[4:
5,]） #10%分位点的收入的消费估计值
qs.rich=c（c（1,x.rich）%*%z$sol[4:
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.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="gray",col.vert="gray"）
ps.wts=（c（0,diff（ps））+c（diff（ps）,0））/2
ap=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="FoodExpenditure",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"））
上图表示收入（income）为10%分位点处（poor，穷人）和90%分位点处（rich，富人）的食品支出的比较。
从左图可以发现，对于穷人而言，在不同分位点估计的食品消费差别不大。
而对于富人而言，在不同分位点对食品消费的差别比较大。
右图反应了穷人和富人的食品消费分布曲线。
穷人的食品消费集中于400左右，比较陡峭；而富人的消费支出集中于800到1200之间，比较分散。
（四）模型比较
# 比较不同分位点下，收入对食品支出的影响机制是否相同
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）
结果：
QuantileRegressionAnalysisofDevianceTable
Model:
foodexp~income
JointTestofEqualityofSlopes:
tauin{ 0.250.50.75 }
DfResidDfFvalue Pr（>F）
1 2 703 15.5572.449e-07***
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
其中P值远小于0.05，故不同分位点下收入对食品支出的影响机制不同。
（五）残差形态的检验
也可以理解为是比较不同分位点的模型之间的关系。
主要有两种模型形式：
（1）位置漂移模型：
不同分位点的估计结果之间的斜率相同或近似，只是截距不同；表现为不同分位点下的拟合曲线是平行的。
（2）位置-尺度漂移模型：
不同分位点的估计结果之间的斜率和截距都不同；表现为不同分位点下的拟合曲线不是平行的。
# 残差形态的检验
source（"C:
/Program"）
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"）
结果：
运行T1，可以查看其检验结果。
其中nullH表示原假设为“location”，即原假设为位置漂移模型。
Tn表示模型整体的检验，统计量为4.8。
THn是对每个自变量的检验。
比较T1和T3的结果（T3的原假设为“位置尺度漂移模型”），T1的统计量大于T3的统计量，可见相对而言，拒绝“位置漂移模型”的概率更大，故相对而言“位置尺度漂移模型”更加合适一些。
>T1
$nullH
[1]"location"
$Tn
[1]4.803762
$THn
X1 X2 X3 X4
1.00031990.53216930.50208340.8926828
attr（,"class"）
[1]"KhmaladzeTest"
>T3
$nullH
[1]"location-scale"
$Tn
[1]2.705583
$THn
X1 X2 X3 X4
1.21028990.69317850.50451630.8957127
attr（,"class"）
[1]"KhmaladzeTest"
（六）非线性分位数回归
这里的非线性函数为Frankcopula函数。
##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# 初始参数
set.seed（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<-data.frame（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
结果：
拟合结果：
（过程略）
>cfMat
分位点 delta mu sigma
[1,] 0.2514.87165-0.205300410.9134657
[2,] 0.5016.25362 0.032325250.9638209
[3,] 0.7512.09836 0.119986140.9423476
（七）半参数和非参数分位数回归
非参数分位数回归在局部多项式的框架下操作起来更加方便。
可以基于以下函数。
#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=0.5）{
# 这是自定义的一个非参数计算函数，在其他数据下同样可以使用
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=0.5） # 关键拟合函数
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=0.75）
fit<-rqss（y~qss（x,lambda=1）,tau=0.9） # 关键的拟合步骤
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=0.5,type="n",main="两组分位数回归结果比较"）
points（engel1,cex=0.5,col=2）
points（engel2,cex=0.5,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（0.05,0.25,0.5,0.75,0.95） # 需要考察的不同分位点
MM2005（foodexp~income,taus,data=engel,group=group,pic=F）# 参数说明，见①
说明：
① 自编函数MM2005的参数说明：
函数调用形式MM2005（formu,taus,data,group,pic=F），其中
#furmu 为回归方程，如foodexp~income；
#taus 为不同的分位数，如taus=c（0.05,0.5,0.95）；
#data 总的数据集，如上例中的engel数据集；
#group 分组指标，是一个向量，用于按行区分data，第一组为1，第二组为2；目前仅能分两组；
#pic 逻辑参数：
是否画图。
如果分位数比较多，建议不画图；默认不画图，pic=F；如果想画图，则添加参数pic=T。
② 最终结果：
>MM2005（foodexp~income,taus,data=engel,group=group,pic=F）
分位数 总差异 回报影响变量影响
0.05-30.452061-72.3593941.90733
0.25 -2.017317-46.2012544.18394
0.50 30.941212-23.2404254.18163
0.75 43.729025-15.7628359.49185
0.95 52.778985-11.2993264.07830
第三节 一个案例
这里使用某年份流动人口中个人消费和个人收入的数据作为案例，进行整体性分析。
（一）初始化程序
基础数据是dta格式（stata软件的数据格式），所以要用foreign包里的read.dta函数读取。
这部分的主要内容包括：
1、修改工作目录至程序和数据所在文件夹，注意不能使用“\”，必须使用“/”或“\\”；
2、清除其他变量、清理内存；
3、加载分位数回归所需的包quantreg和导入数据所需的foreign包；
4、加载4-functions.R文件中的R程序，这里有三个有用的函数；
5、读取案例数据20121218.dta；
6、初始化处理：
去掉含有缺失值的行，检查样本个数，检查列名，整理列名，进行对数变换并生成新的数据集dat2。
# 一个案例
setwd（"F:
/实践2/2012_分位数回归等两个任务/分位数回归方法夹"）# 设置工作目录
setwd（””）
rm（list=ls（））
gc（）
library（quantreg）
library（foreign）
source（"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城镇
# 这里考察上个月的收入与食品支出的关系，并分性别进行讨论
# 把变量的名称改一下
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=na.omit（dat2） # 去掉log变换以后的缺失值行
dim（dat2）
（二）画散点图
画散点
升级会员 