用R语言进行分位数回归

非线性分位数回归这里的非线性函数为Frank copula函数。

 

（六）非线性分位数回归

    这里的非线性函数为Frank copula函数。

## Demo of nonlinear quantile regression model based on Frank copula

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(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(" ")                # 如果前面把每步的结果显示出来，则每次的结果之间添加换行符

}

colnames(cfMat) <- c("分位点",names(coef(fit))) # 给保存系数的矩阵添加列名

cfMat

结果：

图2.5 非参数散点图、真实T分布和拟合结果的比较

拟合结果：（过程略）

> 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

 

（七）半参数和非参数分位数回归

非参数分位数回归在局部多项式的框架下操作起来更加方便。可以基于以下函数。

# 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(i in 1: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(i in hs){

 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.6 基于正态分布为核函数的非参数拟合结果

# 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表示在上图中添加这里的曲线

结果：

图2.7 90%分位点下的附加分位数回归拟合结果

（八）分位数回归的分解

# 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( i in 1: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( i in 1: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.35939 41.90733

   0.25 -2.017317 -46.20125 44.18394

   0.50 30.941212 -23.24042 54.18163

   0.75 43.729025 -15.76283 59.49185

   0.95 52.778985 -11.29932 64.07830