在多元回歸分析中已經介紹過，當自變量之間具有顯著的相關關系時，可能會存在多重共線性。嚴重的多重共線性會大大影響模型的預測結果。除了可以用容忍度與方差擴大因子來度量模型的多重共線性以外，還可以用條件數來度量，常用κ表示，條件數可以定義為：
，
其中，λ為
的特征值（X代表自變量矩陣）。一般認為，當κ>15時，有共線性問題，當κ>30時，說明有嚴重的共線性問題。
本文擬采用R中lars包自帶的diabetes數據集，該數據描述了一些關於糖尿病的血液等化驗指標。通過str()函數先看一下該數據的結構。
> str(diabetes) 'data.frame': 442 obs. of 3 variables: $ x : AsIs [1:442, 1:10] 0.03808 -0.00188 0.0853 -0.08906 0.00538 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : NULL .. ..$ : chr "age" "sex" "bmi" "map" ... $ y : num 151 75 141 206 135 97 138 63 110 310 ... $ x2: AsIs [1:442, 1:64] 0.03808 -0.00188 0.0853 -0.08906 0.00538 ... ..- attr(*, ".Names")= chr "age" "age" "age" "age" ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr "1" "2" "3" "4" ... .. ..$ : chr "age" "sex" "bmi" "map" ...
上說結果說明，diabetes數據為一個數據框，共有3列，442行。第一列是一個442X10數據框，是標准化后的自變量，第二列是因變量，第三列是442X64數據框，包括前者及一些交互作用。
VIF的值可以用軟件包car中的vif()函數計算，條件數κ可以用R內置函數kappa()計算。
> write.csv(diabetes,"diabetes.csv",row.names = FALSE)
> diabetes.test<-read.csv("diabetes.csv")
> kappa(diabetes.test[,12:75])
[1] 11427.09
> vif.dia<-vif(lm(y~.,data = diabetes.test[,11:75]))
> sort(vif.dia,decreasing = TRUE)[1:5]
x2.tc x2.ldl x2.hdl x2.ltg x2.tc.ltg
1295001.21 1000312.11 180836.83 139965.06 61177.87
不管是從κ的值還是從VIF的值來看，該模型都存在嚴重的多重共線性。
1 嶺回歸
假定自變量數據矩陣
為nxp矩陣，通常最小二乘法（ordinary least squares，或ols）通過使殘差平方和達到最小以求得相應的β，即
.
嶺回歸則需要一個懲罰項來約束系數的大小，即在上面的公式中增加
，既要使得殘差平方和小，又要避免系數過大。
說白了，嶺回歸就是通過犧牲最小二乘法的無偏性，以損失部分信息、降低精度為代價獲得回歸系數，對於病態數據以及共線性強的數據具有良好效果。也就是說,在約束條件
下,滿足
.
顯然這里需要確定λ和s的值，一般用交叉驗證或者MallowsCp等准則通過計算來確定。
R中的MASS包提供了lm.ridge()函數用於實現嶺回歸。
> library(MASS) > ridge.dia<-lm.ridge(y~.,lambda = seq(0,150,length.out = 151),data = diabetes.test[,11:75],model = TRUE) #擬合回歸模型 > names(ridge.dia) #查看回歸模型的詳細內容 >[1] "coef" "scales" "Inter" "lambda" "ym" "xm" "GCV" "kHKB" "kLW" >#coef：回歸系數矩陣，每一行對應一個λ。 >#scales：自變量矩陣的標度。 >#Inter：是否包括截距？ >#lambda：λ向量。 >#ym：因變量的均值。 >#xm：自變量矩陣按列計算的均值。 >#GCV：廣義交叉驗證GVC向量。 >#kHKB：嶺常量的kHKB估計。 >#kLW：嶺常量的L-W估計。 > ridge.dia$lambda[which.min(ridge.dia$GCV)] #輸出使GCV最小時的λ值 [1] 81 > ridge.dia$coef[,which.min(ridge.dia$GCV)] #找到GCV最小時對應的系數，該結果是一個向量，長度為自變量的個數 > par(mfrow=c(1,2)) >#繪制回歸系數關於λ的圖形，並作出GCV取取最小值的那條豎線。 > matplot(ridge.dia$lambda,t(ridge.dia$coef),xlab = expression(lambda),ylab = "Coefficients",type = "l",lty = 1:20) > abline(v=ridge.dia$lambda[which.min(ridge.dia$GCV)]) >#繪制GCV關於λ的圖形，並作出GCV取取最小值的那條豎線。 > plot(ridge.dia$lambda,ridge.dia$GCV,type = "l",xlab = expression(lambda),ylab = expression(beta)) > abline(v=ridge.dia$lambda[which.min(ridge.dia$GCV)])
2 lasso回歸
與嶺回歸類似，但是其懲罰項卻不是系數的平方和，而是系數絕對值之和。即在約束條件
下，系數仍然需要滿足最小二乘法的約束條件。
R中的lars包可以用於計算lasso回歸，計算過程如下。
> library(lars) > x=as.matrix(diabetes.test[,1:10]) > y=as.matrix(diabetes.test[,11]) > x2=as.matrix(diabetes.test[,12:75]) > lasso.dia<-lars(x2,y) #lars函數只接受矩陣，所以要把上述的數據框轉化為矩陣 > plot(lasso.dia) #繪制回歸結果，該圖最左邊是只包含截距一種成分，最右邊是包含所有的變量，在不同的步數下，系數變化比較直觀。 > summary(lasso.dia) #輸出lasso對象的細節，包括Df、RSS和Cp值。（Cp是MallowsCp統計量，通常選取Cp最小的那個模型）。
如果從k個自變量中選取p個參與回歸，那么Cp統計量定義為：
其中：
，表示p個自變量參與回歸的誤差平方和。
Ypi：該模型的第i個預測值。
S2：均方殘差，可用樣本的MSE估計。
N：樣本容量。
詳細信息參見：https://en.wikipedia.org/wiki/Mallows%27s_Cp
> cvdia<-cv.lars(x2,y,K = 10) #進行10折交叉驗證
> best<-cvdia$index[which.min(cvdia$cv)] #選擇適合的值（隨機性使得結果不同）
> coef<-coef.lars(lasso.dia,mode="fraction",s=best) #使得CV最小步時的系數。
>該函數同時會繪制CV的變化圖，可以看出在什么時候CV達到極小值。由於在做k折交叉驗證的時候，樣本的分組是隨機地，因此用CV和Cp
>選擇的結果可能會不同
> min(lasso.dia$Cp) #輸出最小的Cp值，結果是第15個
[1] 18.19822
> coef1=coef.lars(lasso.dia,mode="step",s=15) #使lasso.dia最小步時的系數
> CpStep<-data.frame(step=0:104,Df=lasso.dia$df,RSS=lasso.dia$RSS,Cp=lasso.dia$Cp)
> CpStep[9:20,] #可以看出，在第15步時，Cp達到最小值18.2。
step Df RSS Cp
8 8 9 1344199 50.39965
9 9 10 1334100 48.83518
10 10 11 1322126 46.60943
11 11 12 1260433 26.83646
12 12 13 1249726 25.05780
13 13 14 1234993 21.85823
14 14 15 1225552 20.52619
15 15 16 1213289 18.19822
16 16 17 1212253 19.83278
17 17 18 1210149 21.08995
18 18 19 1206003 21.62691
19 19 20 1204605 23.13358
> coef[coef!=0] #ceof!=0,表示變量被選擇，這里根據CV的值選擇了32個變量。
x2.age x2.sex x2.bmi x2.map x2.tc x2.hdl
5.486124 -191.689017 494.870609 301.605908 -47.000781 -248.744630
x2.ltg x2.glu x2.age.2 x2.bmi.2 x2.ltg.2 x2.glu.2
510.705978 50.731562 50.410322 40.870850 -5.528215 93.918944
x2.age.sex x2.age.map x2.age.ldl x2.age.hdl x2.age.ltg x2.age.glu
146.377002 25.700401 -48.088912 26.623122 53.103657 20.494559
x2.sex.bmi x2.sex.map x2.sex.ldl x2.sex.hdl x2.bmi.map x2.map.hdl
31.092360 51.752425 -19.335826 51.477231 124.491470 36.622384
x2.map.glu x2.tc.hdl x2.tc.tch x2.ldl.ltg x2.ldl.glu x2.hdl.tch
-29.354628 8.037888 -60.747576 79.299479 8.188651 -59.171992
x2.tch.ltg x2.tch.glu
-67.607937 24.396984
> coef[coef1!=0] #根據Cp選擇了14個變量
x2.sex x2.bmi x2.map x2.hdl x2.ltg x2.glu x2.age.2
-191.68902 494.87061 301.60591 -248.74463 510.70598 50.73156 50.41032
x2.bmi.2 x2.glu.2 x2.age.sex x2.age.map x2.age.ltg x2.age.glu x2.bmi.map
40.87085 93.91894 146.37700 25.70040 53.10366 20.49456 124.49147
3 適應性lasso回歸
適應性lasso回歸的系數β也要滿足最小二乘公式的條件，但是，其懲罰項為系數絕對值的加權平均值。即約束條件變為：
，
式中，
，而γ是一個調整參數。上面介紹的嶺回歸和lasso回歸是其方法的特例。
R中的msgps包不僅提供了計算適應性lasso（alasso）回歸的方法，還提供了彈性網絡及廣義彈性網絡等方法。其最優策略包括Cp統計量、AIC信息准則、GCV准則以及BIC信息准則。
實現適應性lasso回歸的過程如下。
> #adaptive lasso
> y=diabetes.test[,11]
> al=msgps(x2,y,penalty = "alasso",gamma = 1,lambda = 0)
> summary(al) #
Call: msgps(X = x2, y = y, penalty = "alasso", gamma = 1, lambda = 0)
Penalty: "alasso"
gamma: 1
lambda: 0
df:
tuning df
[1,] 0.000 0.000
[2,] 0.449 1.037
[3,] 1.698 4.412
[4,] 2.957 9.147
[5,] 4.400 12.255
[6,] 5.469 13.958
[7,] 6.677 16.225
[8,] 8.009 18.733
[9,] 9.767 22.160
[10,] 12.113 26.121
[11,] 14.386 28.954
[12,] 16.463 31.118
[13,] 18.832 33.066
[14,] 21.622 36.323
[15,] 24.583 40.062
[16,] 27.704 43.350
[17,] 30.924 45.307
[18,] 34.629 47.100
[19,] 38.619 49.103
[20,] 45.410 54.288
tuning.max: 45.41
ms.coef:
Cp AICC GCV BIC
(Intercept) 1.521e+02 1.521e+02 1.521e+02 152.13
x2.age 0.000e+00 0.000e+00 0.000e+00 0.00
x2.sex -2.138e+02 -2.138e+02 -2.169e+02 -163.49
x2.bmi 5.054e+02 5.054e+02 5.054e+02 505.39
x2.map 3.108e+02 3.108e+02 3.139e+02 275.38
x2.tc -1.486e+02 -1.486e+02 -1.486e+02 -148.57
x2.ldl 1.678e+01 1.678e+01 1.803e+01 10.57
x2.hdl -2.294e+02 -2.294e+02 -2.300e+02 -219.44
x2.tch 0.000e+00 0.000e+00 0.000e+00 0.00
x2.ltg 7.031e+02 7.031e+02 7.031e+02 672.61
x2.glu 0.000e+00 0.000e+00 0.000e+00 0.00
x2.age.2 0.000e+00 0.000e+00 0.000e+00 0.00
x2.bmi.2 0.000e+00 0.000e+00 0.000e+00 0.00
x2.map.2 0.000e+00 0.000e+00 0.000e+00 0.00
x2.tc.2 9.573e+01 9.573e+01 1.392e+02 -36.68
x2.ldl.2 -7.149e+01 -7.149e+01 -1.001e+02 0.00
x2.hdl.2 -2.487e+01 -2.487e+01 -3.170e+01 0.00
x2.tch.2 8.579e+01 8.579e+01 9.387e+01 13.68
x2.ltg.2 3.506e+02 3.506e+02 3.655e+02 159.76
x2.glu.2 3.605e+01 3.605e+01 4.103e+01 0.00
x2.age.sex 1.616e+02 1.616e+02 1.635e+02 99.46
x2.age.bmi 0.000e+00 0.000e+00 0.000e+00 0.00
x2.age.map 0.000e+00 0.000e+00 0.000e+00 0.00
x2.age.tc 0.000e+00 0.000e+00 0.000e+00 0.00
x2.age.ldl 0.000e+00 0.000e+00 0.000e+00 0.00
x2.age.hdl 0.000e+00 0.000e+00 0.000e+00 0.00
x2.age.tch 0.000e+00 0.000e+00 0.000e+00 0.00
x2.age.ltg 3.170e+01 3.170e+01 3.730e+01 0.00
x2.age.glu 0.000e+00 0.000e+00 0.000e+00 0.00
x2.sex.bmi 0.000e+00 0.000e+00 0.000e+00 0.00
x2.sex.map 1.616e+01 1.616e+01 2.487e+01 0.00
x2.sex.tc 2.860e+01 2.860e+01 3.605e+01 0.00
x2.sex.ldl -5.408e+01 -5.408e+01 -6.341e+01 0.00
x2.sex.hdl 2.424e+01 2.424e+01 2.922e+01 0.00
x2.sex.tch 0.000e+00 0.000e+00 0.000e+00 0.00
x2.sex.ltg 0.000e+00 0.000e+00 0.000e+00 0.00
x2.sex.glu 0.000e+00 0.000e+00 0.000e+00 0.00
x2.bmi.map 1.138e+02 1.138e+02 1.181e+02 68.38
x2.bmi.tc 0.000e+00 0.000e+00 0.000e+00 0.00
x2.bmi.ldl 0.000e+00 0.000e+00 0.000e+00 0.00
x2.bmi.hdl 0.000e+00 0.000e+00 0.000e+00 0.00
x2.bmi.tch 0.000e+00 0.000e+00 0.000e+00 0.00
x2.bmi.ltg 0.000e+00 0.000e+00 0.000e+00 0.00
x2.bmi.glu 0.000e+00 0.000e+00 0.000e+00 0.00
x2.map.tc 1.927e+01 1.927e+01 1.927e+01 11.19
x2.map.ldl 0.000e+00 0.000e+00 0.000e+00 0.00
x2.map.hdl 0.000e+00 0.000e+00 7.460e+00 0.00
x2.map.tch 0.000e+00 0.000e+00 0.000e+00 0.00
x2.map.ltg 0.000e+00 0.000e+00 0.000e+00 0.00
x2.map.glu 0.000e+00 0.000e+00 0.000e+00 0.00
x2.tc.ldl 8.889e+01 8.889e+01 8.889e+01 60.92
x2.tc.hdl -3.331e-15 -3.331e-15 -3.331e-15 19.89
x2.tc.tch -2.611e+02 -2.611e+02 -2.741e+02 -135.52
x2.tc.ltg -6.278e+02 -6.278e+02 -6.614e+02 -299.63
x2.tc.glu 0.000e+00 0.000e+00 0.000e+00 0.00
x2.ldl.hdl -9.076e+01 -9.076e+01 -1.113e+02 -23.00
x2.ldl.tch 0.000e+00 0.000e+00 0.000e+00 0.00
x2.ldl.ltg 5.358e+02 5.358e+02 5.446e+02 333.20
x2.ldl.glu 0.000e+00 0.000e+00 0.000e+00 0.00
x2.hdl.tch -2.362e+01 -2.362e+01 -2.362e+01 -23.62
x2.hdl.ltg 2.455e+02 2.455e+02 2.499e+02 128.06
x2.hdl.glu 0.000e+00 0.000e+00 1.243e+00 0.00
x2.tch.ltg 0.000e+00 0.000e+00 0.000e+00 0.00
x2.tch.glu 9.573e+01 9.573e+01 9.573e+01 78.33
x2.ltg.glu 0.000e+00 0.000e+00 0.000e+00 0.00
ms.tuning:
Cp AICC GCV BIC
[1,] 7.987 7.987 8.433 5.111
ms.df:
Cp AICC GCV BIC
[1,] 18.68 18.68 19.47 13.27
>plot(al)
調整參數選的是γ=45.21，為最大值，這時的各個准則及自由度如上所示。
4 偏最小二乘回歸
相關細節可參見：https://en.wikipedia.org/wiki/Partial_least_squares_regression
偏最小二乘回歸有點類似於主成分回歸。偏最小二乘回歸是在因變量和自變量中先各自尋找一個因子，條件是這兩個因子在其他可能的成分中最相關，然后在這選中的一對因子的正交空間中再選一對最相關的因子。如此下去，直到這些對因子有充分代表性為止。R中的pls包可以用來計算偏最小二乘回歸，其計算過程如下。
> ap=plsr(y~x2,64,validation="CV") #求出所以可能的64個因子
> ap$loadings #看代表性，選取一部分因子，部分結果如下圖所示，可以看出，前29個因子的累計方差已經可以解釋總方差的80%以上。
Comp 1 Comp 2 Comp 3 Comp 4 Comp 5 Comp 6 Comp 7 Comp 8 Comp 9
SS loadings 1.142 1.612 1.371 1.378 1.516 1.494 1.476 1.494 1.938
Proportion Var 0.018 0.025 0.021 0.022 0.024 0.023 0.023 0.023 0.030
Cumulative Var 0.018 0.043 0.064 0.086 0.110 0.133 0.156 0.179 0.210
Comp 10 Comp 11 Comp 12 Comp 13 Comp 14 Comp 15 Comp 16 Comp 17
SS loadings 1.890 2.152 1.847 1.693 1.667 1.487 2.532 2.056
Proportion Var 0.030 0.034 0.029 0.026 0.026 0.023 0.040 0.032
Cumulative Var 0.239 0.273 0.302 0.328 0.354 0.377 0.417 0.449
Comp 18 Comp 19 Comp 20 Comp 21 Comp 22 Comp 23 Comp 24 Comp 25
SS loadings 1.854 1.690 1.633 1.785 1.905 2.101 1.579 1.442
Proportion Var 0.029 0.026 0.026 0.028 0.030 0.033 0.025 0.023
Cumulative Var 0.478 0.504 0.530 0.558 0.588 0.620 0.645 0.668
Comp 26 Comp 27 Comp 28 Comp 29 Comp 30 Comp 31 Comp 32 Comp 33
SS loadings 2.058 2.157 1.936 2.010 1.730 1.614 1.865 2.170
Proportion Var 0.032 0.034 0.030 0.031 0.027 0.025 0.029 0.034
Cumulative Var 0.700 0.734 0.764 0.795 0.822 0.847 0.877 0.911
Comp 34 Comp 35 Comp 36 Comp 37 Comp 38 Comp 39 Comp 40 Comp 41
SS loadings 2.286 1.477 1.672 2.145 1.823 1.736 2.173 1.608
Proportion Var 0.036 0.023 0.026 0.034 0.028 0.027 0.034 0.025
Cumulative Var 0.946 0.969 0.995 1.029 1.057 1.085 1.119 1.144
Comp 42 Comp 43 Comp 44 Comp 45 Comp 46 Comp 47 Comp 48 Comp 49
SS loadings 1.664 4.401 1.450 1.796 1.907 2.436 1.267 2.017
Proportion Var 0.026 0.069 0.023 0.028 0.030 0.038 0.020 0.032
Cumulative Var 1.170 1.238 1.261 1.289 1.319 1.357 1.377 1.408
Comp 50 Comp 51 Comp 52 Comp 53 Comp 54 Comp 55 Comp 56 Comp 57
SS loadings 4.171 1.805 1.278 5.221 1.332 10.436 1.160 1.863
Proportion Var 0.065 0.028 0.020 0.082 0.021 0.163 0.018 0.029
Cumulative Var 1.473 1.502 1.522 1.603 1.624 1.787 1.805 1.834
Comp 58 Comp 59 Comp 60 Comp 61 Comp 62 Comp 63 Comp 64
SS loadings 1.302 1.919 4.460 1.134 1.033 1.000 1.000
Proportion Var 0.020 0.030 0.070 0.018 0.016 0.016 0.016
Cumulative Var 1.855 1.885 1.954 1.972 1.988 2.004 2.019
>ap$coefficients #查看各個因子作為原變量的線性組合時的系數，該結果是一個多維數組。
> dim(ap$coefficients)
[1] 64 1 64
> validationplot(ap) #該圖可用來根據CV准則挑選因子數量。可以從圖中看出，5個因子時RMSEP值最小。
致謝：部分內容來自於吳喜之《復雜數據分析——基於R的應用》，數據為R中lars包自帶數據集diabetes。