用R语言进行分位数回归Word文档格式.docx
《用R语言进行分位数回归Word文档格式.docx》由会员分享，可在线阅读，更多相关《用R语言进行分位数回归Word文档格式.docx（35页珍藏版）》请在冰豆网上搜索。
（4）
样本分位数回归是使加权误差绝对值之和最小，即
（5）
上式可等价表示为：
其中，
为检查函数（checkfunction），定义为：
为指示函数（indicatorfunction），z是条件关系式，当z为真时，
；
当z为假时，
。
同线性方程y=kx比较，
相当于直线的斜率k，可以看出，
为分段函数，如下图所示。
现假设因变量Y由k个自变量组成的矩阵X线性表示，对于条件均值函数
，通过求解（8）式得到参数估计值
对于条件分位数函数，通过求解（9）式得到参数估计值
式中，
函数表示取函数最小值时
的取值。
（三）分位数回归模型的参数估计算法
1、主要算法
（1）单纯形算法（SimplexMethod）
Koenker和Orey[5]（1993）把分两步解决最优化问题的单纯形算法[6]扩展到所有回归分位数中。
该算法估计出来的参数具有很好的稳定性，但是在处理大型数据时运算的速度会显著降低。
（2）内点算法（InteriorPointMethod）
由于单纯形算法在处理大型数据时效率低下，Karmarker提出了内点算法[7]；
Portnoy和Koenker把这种方法是用在分位数回归中，得出了处理大型数据时内点算法的运算速度远快于单纯形算法的结论。
但内点算法每计算一步都要进行因数分解，当自变量比较多的时候效率比较低。
其次，如果要达到和单纯形算法一样的精度，就必须进行舍入步骤的计算，者也降低了算法的运行效率。
（3）平滑算法（SmoothingMethod）
上述两种算法都有各自的优点和不足，而有限平滑算法则是一种同时兼顾运算效率以及运算速度的方法。
Chen把这种算法扩展到计算回归分位数中[8]。
2、R语言quantreg包中的假设检验
加载quantreg包以后，使用summary（）函数或summary.rq（）函数，可以得到参数系数的一些假设检验统计量。
其实，以上两个函数是一致的。
在使用summary（）的时候，如果sumamry（）加载的模型（对象）是分位数回归模型，则会自动调用summary.rq（）来处理这个对象。
summary.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不同，这里是参数se的一个取值，而上面的iid是一个逻辑参数）假定残差服从独立同分布，并按照KB（1978）的方法计算残差。
nid:
用sparsity算法计算的参数估计值标准差。
ker:
用Powell（1990）的核密度估计方法得到标准差。
boot:
采用bootstrap自助抽样的方法计算标准差。
默认情况下，se=NULL且convariance=FALSE，标准差的默认算法是se=”rank”；
其他情况下，se默认值为”nid”。
#covariance:
逻辑参数，是否返回参数估计量的协方差矩阵。
不同参数的结果，可参看下面的程序案例。
（四）分位数分解（MM2005方法）[9]
我们可以进一步运用分位数分解法对各个影响因素进行分解分析[10]。
这里仅介绍MM2005方法。
为讲解方便，这里以各因素对城乡家庭收入的影响为例，观察各个影响因素在不同分位数上对城乡家庭收入差异的影响度的大小。
这里介绍Machado和Mata[11]（2005）提出的分位数分解法，将每个分位数上的城乡收入差异分解为两个部分：
一部分是由于城乡家庭劳动力特征的不同回报率引起的（即分位数回归参数的不同引起的，TheReturnEffects），例如城乡家庭劳动力在相同的教育程度、工作年限以及所处当地的经济发展水平相同的特定因素下不同的回报率引起的家庭人均收入差异；
另一部分是由于城乡家庭劳动力的特征变量分布不同引起的（即影响因素变量值的不同引起的，TheCovariateEffect），城乡家庭人均收入这部分的差异会随着样本分布的不同而略有变化。
利用Machado和Mata分位数分解方法的关键是进行反事实分析（thecounter-factualanalysis），我们最关心的一种反事实分析就是，如果城市家庭劳动力按照农村家庭劳动力的分位数回归参数决定家庭人均收入的话，城市家庭的人均收入分布会如何？
这里定义反事实分布为
，其中
表示影响城市家庭人均收入的变量分布，
表示影响农村家庭人均收入的变量在每个分位数上的回归参数。
表示如果城市家庭劳动力按照农村家庭劳动力的分位数回归参数决定家庭人均收入的话，城市家庭的反事实人均收入的大小。
的具体计算步骤为：
（1）确定不同的分位点，分别表示为
（2）在农村家庭样本中，分别以
做分位数回归，得到
组分位数回归参数向量
（3）将城市家庭样本数据表示为
（4）把
（2）中得到的分位数回归参数和（3）中得到得城市家庭子样本变量分布相结合，得到一个新的样本，即反事实分布样本
假定在τ分位数下城市家庭人均收入、反事实家庭人均收入和农村家庭人均收入分别为
、
则不同分位数下的城乡家庭人均收入分布差异可表示为：
等式右边的第一项称为“回报影响（thereturneffect）”，它表示在不同的分位数下，由于城乡家庭劳动力的生产回报率不同所导致的城乡差异部分；
等式右边第二项成为“变量影响（thecovariateeffect）”，它表示不同分位数下城乡家庭随机抽样的样本变量分布不同所导致的城乡差异部分。
（五）非线性分位数回归和非参数分位数回归
暂略。
第二节
（一）安装和加载包
R语言的基本包中没有进行分位数回归的程序包，故需要在官网下载并安装相应的程序包quantreg。
在电脑上安装过quantreg包以后，下次不需要再次安装了。
但每次使用分位数回归前，需要加载quantreg包。
install.package（“quantreg”）
#保持联网的情况下安装包
library（“quantreg”）
#加载包
help.start（）
#进入R帮助首页
help（rq）
#获取rq函数的帮助，也可以写成：
?
rq
example（rq）
#显示分位数回归函数rq（）的一个简单示例代码
（二）一个简单的分位数回归模型及结果
data（engel）
#加载quantreg包自带的数据集，见说明①
fit1=rq（foodexp~income,tau=0.5,data=engel）
#进行分位数回归，见说明②
fit1
#直接显示分位数回归的模型和系数，见说明③
summary（fit1）
#得到更加详细的显示结果，见说明④
r1=resid（fit1）
#得到残差序列，并赋值为变量r1
c1=coef（fit1）
#得到模型的系数，并赋值给变量c1，见说明⑤
summary（fit1,se=“nid”）
#通过设置参数se，可以得到系数的假设检验，说明⑥
说明：
①engel（1857）是考察食物支出与家庭收入之间关系的一个数据集，用函数head（engel）可以查看前六行的值：
#
income
foodexp
#1420.1577255.8394
#2541.4117310.9587
#3901.1575485.6800
#4639.0802402.9974
#5750.8756495.5608
#6945.7989633.7978
②这里因变量为foodexp，即食物支出。
自变量为income，即家庭收入。
-tau表示计算50%分位点的参数，这里可以同时计算多个分位点的分位数回归结果，如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”，针对特别大数据，使用经过预处理的Frisch–Newton逼近方法；
D.“fnc”，针对被拟合系数特殊的线性不等式约束情况；
E.“lasso”和“scad”，基于特定惩罚函数的平滑算法进行拟合。
③直接运行fit1，会得到简单的计算结果，如：
#Call:
rq（formula=foodexp~income,tau=0.5,data=engel）
#
#Coefficients:
（Intercept）
income
81.4822474
0.5601806
#Degreesoffreedom:
235total;
233residual
④用summary（）函数可以得到回归模型的详细结果，包括系数和上下限。
#tau:
[1]0.5
coefficientslowerbd
upperbd
#（Intercept）
81.48225
53.25915114.01156
#income
0.56018
0.48702
0.60199
⑤coef（）函数得到的系数为向量形式，第一个元素为常数项的系数，第二个及以后为自变量的系数。
⑥summary函数se参数的说明：
A.se=“rank”:
按照Koenker（1994）的排秩方法计算得到的置信区间，默认残差为独立同分布。
注意的是，上下限是不对称的。
Coefficients:
（Intercept）
income
B.se=”iid”:
假设残差为独立同分布，用KB（1978）的方法计算得到近似的协方差矩阵。
Value
Std.Errortvalue
Pr（>
|t|）
（Intercept）81.4822513.23908
6.15468
0.00000
0.01192
46.99766
C.se=“nid”:
表示按照Huber方法逼近得到的估计量。
（Intercept）81.4822519.25066
4.23270
0.00003
0.02828
19.81032
D.se=”ker”:
采用Powell（1990）的核估计方法。
（Intercept）81.4822530.21532
2.69672
0.00751
0.03732
15.01139
E.se=”boot”:
采用bootstrap方法自助抽样的方法估计系数的误差标准差。
（Intercept）81.4822525.23647
3.22875
0.00142
0.03194
17.53752
（三）不同分位点下的回归结果比较
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）
plot（fit2）
结果：
图2.1
99个分位点的系数估计值
图2.2
5个分位点的系数估计值
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="
）
#画中位数回归的拟合直线，颜色蓝
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）
图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$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:
#90%分位点的收入的消费估计值
windows（10,5）
par（mfrow=c（1,2））
#把绘图区域划分为一行两列
plot（c（ps,ps）,c（qs.poor,qs.rich）,type="
#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="
col.vert="
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="
xlab="
Density"
lines（qs.rich,ar$dens,col="
lines（qs.poor,ap$dens,col="
black"
legend（"
topright"
c（"
poor"
"
rich"
）,lty=c（1,1）,
col=c（"
））
图2.4
10%分位点和90%分位点之间的比较
上图表示收入（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:
/ProgramFiles/R/R-2.15.0/library/quantreg/doc/gasprice.R"
x=gasprice
n=length（x）
p=5
X=cbind（x
升级会员 