我读硕士老师给我的第一篇论文就是一个分位数回归的文章,当时觉得这个模型很简单,我很快就用R的示例文件写了一个例子,但是,在后面的研究中,我越来越觉得,这个模型没有我想的那么简单,而且有着非常丰富的内涵需要来挖掘,就找了好几本书来看,结果真的是越看越懵,越看越懵,但是懵了一段时间之后,又重新感觉自己明白点了,所以赶紧把这一点进行一个总结,省的再放一段时间,连仅有的这一点懂的东西都没有了。
首先随机变量X的分布函数为:
则它的 分位数可定义为:
若将分布函数F(x)的逆定义为:
则
设有随机向量(X,Y),其中Y在给定X=x的情况下的条件累积分布函数为 ,则将该条件随机变量Y|X=x的
条件分位数定义为:
假定我们有样本序列 满足下列回归模型,即
假定误差项 为独立同分布的序列,且分布情况未知,这Y的
条件分位数
式中, 为分位回归领域的损失函数,有些地方也称之为检验函数。
为示性函数。不包含示性函数的损失函数为
损失函数同样可以用下面一个公式进行表示:
则我们希望 尽可能的小,最好等于0,则:
最优解为
用样本经验分布函数代替F(x)则
这些是分位数回归的一些基本思想,这仅仅只涉及到基本的分位数线性回归,还有较为复杂的核回归在随后的研究中再进行更新。
下面我们来看一下代码实现。
分位数回归可以使用R中的quantreg包进行实现,我们就是用该包中的示例数据构建分位数回归模型:
首先加载数据和包
library(SparseM)
library(quantreg)
data("engel")
engel的这个数据中有两个变量一个变量是income,代表每个家庭的收入,另一个变量是foodexp代表食物支出,数据的大致情况是这样的:
![](https://i-blog.csdnimg.cn/blog_migrate/ec84eace3abd9d7455852c8836f4920f.png)
我们来建立一个0.5分位数回归:
fit1 <- rq(foodexp ~ income, tau = .5, data = engel)
fit1
summary(fit1)
我们使用summary来看一下模型的大致情况:
![](https://i-blog.csdnimg.cn/blog_migrate/e6b2c5b3beadfe5077fb12b6afee0274.png)
大体上来看,income的 值是有意义的(95%置信区间未包含0),即income变动一个单位,foodexp的0.5分位数变动0.56个单位。我们可以通过summary的se参数数据的不同假设,iid即为独立同分布的假设。
我们可以使用quantreg内置的函数计算模型拟合的残差,并绘制残差图:
r1<-resid(fit1)#计算残差
plot(r1)
![](https://i-blog.csdnimg.cn/blog_migrate/3b8133500e3d8f5706fbe10f59fab016.png)
残差图数据均匀分布在(-200,200)的区间中,模型拟合的结果还是可以的。
接下来,我们可以画出,income和foodexp的散点图,并绘制不同分位数的分位回归曲线。
attach(engel)
plot(income, foodexp, cex = 0.25, type = "n", xlab ="Household Income", ylab = "Food Expenditure")
points(income,foodexp , cex = 0.25, col = "blue")
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 (i in 1:length(taus)) {
abline(rq(foodexp ~income , tau = taus[i]), col = "gray")
}
detach(engel)
![](https://i-blog.csdnimg.cn/blog_migrate/4b1eb336332b4d562bf6d531c53bd746.png)
图中的蓝色实线代表0.5分位回归的回归曲线,红色虚线代表使用最小二乘法拟合的线性回归曲线,几条灰色的曲线,代表0.05, 0.1, 0.25, 0.75, 0.9, 0.95的分位数回归曲线。
我们可以看到不同分位数的分位回归曲线的斜率不同,0.5分位数的分位回归曲线和线性回归曲线的斜率也不相同,则可以认为,income对于foodexp不同的分位数影响不同,据此,我们将绘制不同分位数的分位回归曲线的置信区间图,来进一步探讨income对于foodexp的影响。
engel<-within(engel,xx <- income - mean(income))
fit1 <- summary(rq(foodexp~xx,tau=2:98/100,data = engel))
plot(fit1,mfrow = c(1:2))
![](https://i-blog.csdnimg.cn/blog_migrate/04d1c08a46fb971ef2b3810c40cdfc32.png)
上图绘制了0.02到0.98这个区间中,每隔0.01做一次分位回归,其中黑色实心点代表回归曲线的 值,阴影部分代表95%置信区间,红色实线和虚线分别代表的是,线性回归曲线的
值和置信区间。从图中可以看出,income对于0.02分位的foodexp的影响在0.35左右,对于0.98分位的影响在0.7左右。
不同分位数的分位回归的 值是否是由于抽样误差造成的,我们同样需要假设检验进行验证,那么我们使用Wald 检验进行验证。 得到如下结果,P值小于0.05可以认为不同分位数回归的
值是有差异的。
fit1 <- rq(y ~ x, tau = 0.25)
fit2 <- rq(y ~ x, tau = 0.5)
fit3 <- rq(y ~ x, tau = 0.75)
anova(fit1, fit2, fit3)
![](https://i-blog.csdnimg.cn/blog_migrate/7caa73a3f6118bd57f7baea147143043.png)