用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("\n")                # 如果前面把每步的結果顯示出來，則每次的結果之間添加換行符

}

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