用R語言進行分位數迴歸

阿新 • • 發佈：2018-12-04

非線性分位數迴歸這裡的非線性函式為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

用R語言進行分位數迴歸

用R語言進行分位數迴歸

用R語言進行分位數迴歸：基礎篇

用R語言進行分位數回歸：基礎篇

用R語言進行迴歸分析

用R語言進行方差分析

用R語言進行基本統計分析

乾貨：用R語言進行資料提取的方法！

用R語言進行關聯分析

用R語言的quantreg包進行分位數回歸

使用R語言進行中文分詞

用R語言對網路資料進行統計分析(五)

（轉）python中調用R語言通過rpy2 進行交互安裝配置詳解

R語言進行機器學習方法及實例

用R語言分析與預測員工離職

R語言學習筆記(二)——迴歸分析

R語言開發之線性迴歸瞭解下

使用R語言進行時間序列（arima，指數平滑）分析

【機器學習】R語言進行機器學習方法及例項

用R語言實現對不平衡資料的四種處理方法

分位數迴歸-Quantile regression

用R語言進行分位數迴歸

相關推薦