分位数回归.

2、不同分位点拟合曲线的比较# 散点图attach(engel) # 打开engel数据集，直接运行其中的列名，就可以调用相应列plot(income,foodexp,cex=0.25,type="n", # 画图，说明①xlab="Household Income", ylab="Food Expenditure")points(income,foodexp,cex=0.5,col="blue") # 添加点，点的大小为0.5abline( 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)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) )/2ap = 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="Food Expenditure", 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结果：Quantile Regression Analysis of Deviance TableModel: foodexp ~ incomeJoint Test of Equality of Slopes: tau in { 0.25 0.5 0.75 }Df Resid Df F value Pr(>F)1 2 703 15.557 2.449e-07 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘’ 1其中P值远小于0.05，故不同分位点下收入对食品支出的影响机制不同。

（五）残差形态的检验也可以理解为是比较不同分位点的模型之间的关系。

主要有两种模型形式：（1）位置漂移模型：不同分位点的估计结果之间的斜率相同或近似，只是截距不同；表现为不同分位点下的拟合曲线是平行的。

（2）位置-尺度漂移模型：不同分位点的估计结果之间的斜率和截距都不同；表现为不同分位点下的拟合曲线不是平行的。

# 残差形态的检验source("C:/Program Files/R/R-2.15.0/library/quantreg/doc/gasprice.R")x = gaspricen = length(x)p = 5X = 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")结果：运行T1，可以查看其检验结果。

其中nullH表示原假设为“location”，即原假设为位置漂移模型。

Tn表示模型整体的检验，统计量为4.8。

THn是对每个自变量的检验。

比较T1和T3的结果（T3的原假设为“位置尺度漂移模型”），T1的统计量大于T3的统计量，可见相对而言，拒绝“位置漂移模型”的概率更大，故相对而言“位置尺度漂移模型”更加合适一些。

> T1$nullH[1] "location"$Tn[1] 4.803762$THnX1 X2 X3 X41.0003199 0.5321693 0.5020834 0.8926828attr(,"class")[1] "KhmaladzeTest"> T3$nullH[1] "location-scale"$Tn[1] 2.705583$THnX1 X2 X3 X41.2102899 0.6931785 0.5045163 0.8957127attr(,"class")[1] "KhmaladzeTest"（六）非线性分位数回归这里的非线性函数为Frank copula函数。

## Demo of nonlinear quantile regression model based on Frank copulavFrank <- 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(i in 1: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(i in 1: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.25 14.87165 -0.20530041 0.9134657[2,] 0.50 16.25362 0.03232525 0.9638209[3,] 0.75 12.09836 0.11998614 0.9423476（七）半参数和非参数分位数回归非参数分位数回归在局部多项式的框架下操作起来更加方便。