分位数回归（quantile regression）

分位数回归模型的参数估计，除了频率学派的统计方法，如内点法、极大似然估计法之外，还可以运用贝叶斯学派的参数估计法。与频率统计方法相比，贝叶斯参数估计法具有一些独特的优势，在一定程度上弥补了频率统计方法的不足：一是在小样本数据下，频率统计方法面临着由于方差过大而使估计量的优良统计性质难以保障的困难，进而使得以渐进正态性为基础的假设检验失效，但是贝叶斯参数估计法可以避免对渐进协方差的估计，从而保证参数估计值的性质稳定。二是频率统计方法只能得到一个参数估计值。但贝叶斯估计法把参数视作一个随机变量，考虑了参数的不确定风险，提高了估计和预测的精度。

用R语言实现。

因变量Y不服从正态分布

library(quantreg)

data(engel)

plot(engel$income, engel$foodexp, xlab='income', ylab='foodexp')

boxplot(engel$foodexp, xlab='foodexp')

qqnorm(engel$foodexp, main='QQ plot')

qqline(engel$foodexp, col='red', lwd=2)

分位数回归

模型拟合效果

attach(engel)

fit1 <- rq(foodexp ~ income, tau = .5, data = engel)

fit1

summary(fit1)

r1<-resid(fit1)#计算残差

plot(r1)

残差图数据均匀分布在（-200,200）的区间中，模型拟合的结果还是可以的。

绘制不同分位数的分位回归曲线

rq_result <- rq(foodexp ~ income, tau=c(0.05, 0.25, 0.5, 0.75, 0.95))

plot(income, foodexp, cex=0.25, type='n', xlab='income', ylab='foodexp')

points(income, foodexp, cex=0.5, 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')

}

图中的蓝色实线代表0.5分位回归的回归曲线，红色虚线代表使用最小二乘法拟合的线性回归曲线，几条灰色的曲线，代表0.05, 0.1, 0.25, 0.75, 0.9, 0.95的分位数回归曲线。

绘制不同系数的不同分位数置信区间图

engel<-within(engel,xx <- income - mean(income))

fit1 <- summary(rq(foodexp~xx,tau=2:98/100,data = engel))

plot(fit1,mfrow = c(1:2))

上图绘制了0.02到0.98这个区间中，每隔0.01做一次分位回归，其中黑色实心点代表回归曲线的 β 值，阴影部分代表95%置信区间，红色实线和虚线分别代表的是，线性回归曲线的 β值和置信区间。

不同分位数的分位回归的 β值是否是由于抽样误差造成的，我们同样需要假设检验进行验证，那么我们使用Wald 检验进行验证。

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)

面板分位数回归

将截面分位数回归模型扩展至面板数据框架下，考虑了研究单元之间的个体异质性。

library(rqpd)

library(plm)

data(Produc,package="plm")

mydata = Produc

m <- 17    ##17年

n <- 48    ##美国48州

s <- as.factor(rep(1:n,rep(m,n)))   ##声明面板数据结构变量

x <-cbind(mydata$pcap,mydata$hwy)   ##数据中取出解释变量x

y <- mydata$unemp                   ##数据中取出被解释变量y

##rqpd直接运算公式，个体效应惩罚系数lambda为1,系统默认分位点是0.25,0.5,0.75

fit <- rqpd(y ~x | s, panel(lambda = 1))

sfit <- summary(fit)

sfit

sfit$coefficients  ##能一并取出估计的个体效应系数

##依次定义九个分位点同时定义tau和tauw

fit1 <- rqpd(y ~x | s, panel(taus=seq(0.1,0.9,by=0.1),tauw=rep(1/9,9),lambda = 1))

sfit1 <- summary(fit1)

sfit1