在R內實踐嶺回歸 (Doing Ridge regression in R Statistical software)

    如上次的百分位數回歸的一節所說, 作線性回歸的基本要求是LINE, 其中I(數據獨立性)和N(數據呈正態分佈)是較難符合要求的...
    不符合N時, 可試用百分位數回歸方法解決; 而不符合I時, 可試用嶺(岭)回歸方法!
    為方便驗證是操作的正確性, 現引用高惠璇老師在《實用統計方法與SAS系統》P100-103例子作為對照! 好! 現在開始...
有這樣的一個數據:


#在R軟件內輸入並合成數據集d452, x1:國內生產總值, x2:儲存量, x3:總消費量, y:進口總額
x1<-c(149.3,161.2,171.5,175.5,180.8,190.7,202.1,212.4,226.1,231.9,239.0)
x2<-c(4.2,4.1,3.1,3.1,1.1,2.2,2.1,5.6,5.0,5.1,0.7)
x3<-c(108.1,114.8,123.2,126.9,132.1,137.7,146.0,154.1,162.3,164.3,167.6)
y<-c(15.9,16.4,19.0,19.1,18.8,20.4,22.7,26.5,28.1,27.6,26.3)
d452<-data.frame(x1,x2,x3,y)
lm.1<-lm(y~.,data = d452) #試作普通的線性回歸分析
summary(lm.1)
#結果

Call: lm(formula = y ~ ., data = d452)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.52367 -0.38953  0.05424  0.22644  0.78313 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -10.12799    1.21216  -8.355  6.9e-05 ***
x1           -0.05140    0.07028  -0.731 0.488344    
x2            0.58695    0.09462   6.203 0.000444 ***
x3            0.28685    0.10221   2.807 0.026277 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4889 on 7 degrees of freedom
Multiple R-squared:  0.9919, Adjusted R-squared:  0.9884 
F-statistic: 285.6 on 3 and 7 DF,  p-value: 1.112e-07

x1為負值,與現實世界不符合... 

library(car)
vif.1<-vif(lm.1)
vif.1
#結果

        x1         x2         x3 
185.997470   1.018909 186.110015 

x1及x3的VIF大於10,有多重共線性 

cor.1<-cor.test(x1,x3)
cor.1
#結果

 Pearson's product-moment correlation

data:  x1 and x3
t = 40.448, df = 9, p-value = 1.718e-11
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9890918 0.9993142
sample estimates:
      cor 
0.9972607

x1與x3的相關系數達0.997(p<0.01),有強烈的相關性 

#在這種情況下,線性回歸的結果受多重共線性影響,方差會變得很大(因合力問題)!

為使方差變小,得加入一個系數(lambda).但這個系數的最小取值得依賴於方程中的β和σ2.

系數的估計方法有:

1. 嶺跡圖

2. 方差膨脹因子法

3. 控制殘差平方和法

library(MASS) #嶺回歸法在MASS的統計包內---lm.ridge命令
lm.r1<-lm.ridge(y~.,data = d452,lambda = seq(0,0.5,0.01)) #lambda系數可用seq序列, 讓函數自已試
#結果

                          x1        x2        x3
0.00 -10.127988 -0.051396160 0.5869490 0.2868487
0.01  -9.861266 -0.020071030 0.5920564 0.2411973
0.02  -9.696169 -0.001390176 0.5948969 0.2139345
...
0.50  -8.466721  0.064243563 0.5832040 0.1140136

嶺跡圖, 與lm.ridge的列表對應

plot(lm.r1)
select(lm.r1)
#結果

modified HKB estimator is 0.006855363
modified L-W estimator is 0.01283802 
smallest value of GCV  at 0.01

被選出最小的lambda系數是0.01 

lm.r2<-lm.ridge(y~.,data = d452,lambda = 0.01)
lm.r2
#結果

                     x1          x2          x3 
-9.86126637 -0.02007103  0.59205639  0.24119727 

若以lambda系數為0.01再代入計算,則得嶺回歸方程為上.

lm.r3<-lm.ridge(y~.,data = d452,lambda = 0.02)
lm.r3#結果

                       x1           x2           x3 
-9.696168972 -0.001390176  0.594896865  0.213934536

若以該書的0.02代入嶺回歸公式,得上述結果!與書的結果不符合...

lm.r4<-lm.ridge(y~.,data = d452,lambda = 0.22)
lm.r4 

#結果

                     x1          x2          x3 
-8.92769761  0.05734309  0.59541713  0.12663337 

這才與書上的結果一致... 

在R內實踐百分位數回歸 (Doing Quantile Regression in R statistical software)

    近日不約而同地有3位碩士研究生, 都問了類似的問題: “老師, 我想做回歸分析, 而應變項雖然是連續性變項, 但郤是偏態的, 怎麼辦好?”

    於是我都要她們複習一下線性回歸的基本要求, 握要地說是”LINE”, 要求變項:

1. “L”: 是線性的 (連續性數據)
2. “I”: 變項間是互相獨立 (不要有多重共線性)
3. “N”: 要正態分佈
4. “E”: 誤差獨立分佈

    現實中”L”和”E” 基本都是滿足的, 但”I”和”N”就未必了… 人們是睿智的, 想到解決”I”和”N”的問題是:

1. “I”---岭回歸, 主成分分析, 或只選擇最相關的變項作回歸分析;
2. “N”---百分位數回歸…

現在主要介紹百分位數回歸! 其本質是求解任何分位值時, 各引數 x 的回歸係數。其能在一定程度上全面展示所有資料的資訊，但更側重於對特徵分佈不同位置人群的異質性分析。它不僅可以度量回歸變數在分佈中心的影響，而且還可以度量在分佈上尾和下尾的影響，因此較之傳統線性回歸模型具有獨特的優勢，能夠得到全面的分析結果，尤其適用於不均勻的條件分佈，且分位數回歸係數較傳統線性回歸模型的係數估計更穩健。

進行百分位數回歸的統計軟件主要有SAS, R, STATA等, 但以SAS及R為普遍! 在SAS內以proc quantreg為操作步驟. 在R內則以quantreg統計包為主. 現以R為介紹, 並使用quantreg統計包內自帶的數據engel示範…

library(quantreg)  #安裝及使用quantreg統計包 : install.packages(“quantreg”)

data("engel")  #使用自帶的數據engel

food.nor<-shapiro.test(engel$foodexp)  #數據庫內foodexp的正態性檢驗

food.hist<-hist(engel$foodexp)  #作個直條圖

   Shapiro-Wilk normality test

data:  engel$foodexp

W = 0.87471, p-value = 5.45e-13

(明顯地是偏態分佈啦...)

fit=rq(foodexp~income,data=engel,tau=c(0.25,0.5,0.75))  #百分位數回歸: 命令rq, 公式foodexp~income, 使用數據庫engel, 百分位數0.25下四分位數,0.5中位數,0.75上四分位數

summary(fit)  #算出總結果

plot(summary(fit))  #對結果作圖

 

attach(engel)  #打開engel資料集，直接運行其中的列名，就可以調用相應列

plot(income,foodexp,cex=0.25,type="n", xlab="家庭收入", ylab="食物開支")  # 畫圖

points(income,foodexp,cex=0.5,col="blue")  # 添加點，點的大小為0.5

abline( rq(foodexp ~ income, tau=0.5), col="blue" )  # 畫中位數回歸的擬合直線，顏色藍

abline( lm(foodexp ~ income), lty = 2, col="red" )  # 畫普通最小二乘法擬合直線，顏色紅

taus<-c(0.25, 0.5, 0.75)

for(i in 1:length(taus)){  # 繪製不同分位點下的擬合直線，顏色為灰色

  abline( rq(foodexp ~ income, tau=taus[i]), col="gray" )

}

detach(engel)

 參考資料:

研究健康相关生命质量影响因素的分位数回归分析. 中国卫生统计2016年4月第33卷第2期.P190-193.

分位数回归及应用简介. 统计与信息论坛. 第21卷第3期2006年5月.P35-38.