###分位數(shù)回歸####

library(quantreg)??# 載入quantreg包

data(engel)? ?? ???# 加載quantreg包自帶的數(shù)據(jù)集



##單變量回歸

#分位數(shù)回歸(tau = 0.5)

fit1 = rq(foodexp ~ income, tau = 0.5, data = engel)? ?? ?? ?

r1 = resid(fit1)? ?# 得到殘差序列，并賦值為變量r1

c1 = coef(fit1)? ? # 得到模型的系數(shù)，并賦值給變量c1

summary(fit1)? ?? ?# 顯示分位數(shù)回歸的模型和系數(shù)

summary(fit1, se = "boot") # 通過設置參數(shù)se，可以得到系數(shù)的假設檢驗



#分位數(shù)回歸(tau = 0.5、0.75)

fit1 = rq(foodexp ~ income, tau = 0.5, data = engel) 

fit2 = rq(foodexp ~ income, tau = 0.75, data = engel)

#模型比較

anova(fit1,fit2)? ? #方差分析


? ???
? ?? ? col=c("#EEEE00","#EE6363"))

#不同分位點的回歸比較

fit = rq(foodexp ~ income,??tau = c(0.05,0.25,0.5,0.75,0.95), data = engel)

plot( summary(fit))

##多變量回歸

data(barro)

fit1 <- rq(y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro,tau=.25)

fit2 <- rq(y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro,tau=.50)

fit3 <- rq(y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro,tau=.75)

# 替代方式 fit <- rq(y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2, method = "fn", tau = 1:4/5, data = barro)

anova(fit1,fit2,fit3)? ?? ?? ?? ? # 不同分位點模型比較-方差分析

anova(fit1,fit2,fit3,joint=FALSE)

#畫圖比較分析

plot(engelfooder,engelincome,pch=20, col = "#2E8B57",

? ???main = "家庭收入與食品支出的分位數(shù)回歸",xlab="食品支出",ylab="家庭收入")

lines(fitted(fit1), engel$income,lwd=2, col = "#EEEE00")

lines(fitted(fit2), engel$income,lwd=2, col = "#EE6363")

legend("topright", c("tau=.5","tau=.75"), lty=c(1,1),

? ?? ? col=c("#EEEE00","#EE6363"))
#不同分位點的回歸比較

fit = rq(foodexp ~ income,??tau = c(0.05,0.25,0.5,0.75,0.95), data = engel)

plot( summary(fit))

##多變量回歸

data(barro)

fit1 <- rq(y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro,tau=.25)

fit2 <- rq(y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro,tau=.50)

fit3 <- rq(y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2, data = barro,tau=.75)

# 替代方式 fit <- rq(y.net ~ lgdp2 + fse2 + gedy2 + Iy2 + gcony2, method = "fn", tau = 1:4/5, data = barro)

anova(fit1,fit2,fit3)? ?? ?? ?? ? # 不同分位點模型比較-方差分析

anova(fit1,fit2,fit3,joint=FALSE)

#不同分位點擬合曲線的比較

plot(barro$y.net,pch=20, col = "#2E8B57",

? ???main = "不同分位點擬合曲線的比較")

lines(fitted(fit1),lwd=2, col = "#FF00FF")

lines(fitted(fit2),lwd=2, col = "#EEEE00")

lines(fitted(fit3),lwd=2, col = "#EE6363")

legend("topright", c("tau=.25","tau=.50","tau=.75"), lty=c(1,1),

? ?? ? col=c( "#FF00FF","#EEEE00","#EE6363"))