数据分析中常见的七种回归分析以及R语言实现（三）---岭回归

我们在回归分析的时候，古典模型中有一个基本的假定就是自变量之间是不相关的，但是如果我们在拟合出来的回归模型出现了自变量之间高度相关的话，可能对结果又产生影响，我们称这个问题为多重共线性，多重共线性又分为两种，一种是完全多重共线性，还有一种是不完全多重共线性；

产生的原因有几个方面

1、变量之间存在内部的联系

2、变量之间存在共同的趋势等

造成的后果分两部分

完全多重共线性造成的后果

1、当自变量线性相关的时候，参数将无法唯一确定，参数的方差将趋近于无穷大，这时候无法使用最小二乘法

不完全多重共线性造成的后果

1、参数估计量的方差随着多重共线性的严重程度的增加而增加，但是参数是可以估计的

2、进行统计检验时容易删除掉重要解释变量

因为当多重共线性的时候容易造成自变量对因变量不显著，从模型中错误的剔除，这样容易删除重要解释变量的设定；

3、参数的置信区间明显扩大

因为由于存在多重共线性。我们的参数估计都有较大的标准差，因此参数真值的置信区间也将增大

那么我们怎么去判断一个模型上存在多重共线性呢？

根据经验表明，多重共线性存在的一个标志就是就模型存在较大的标准差，和较小的T统计量，如果一个模型的可决系数R^2很大，F检验高度限制，但偏回归系数的T检验几乎都不显著，那么模型很可能是存在多重共线性了。因为通过检验，虽然各个解释变量对因变量的共同影响高度显著，但每个解释变量的单独影响都不显著，我们无法判断哪个解释变量对被解释变量的影响更大

1、可以利用自变量之间的简单相关系数检验

这个方法是一个简便的方法，一般而言，如果每两个解释变量的简单相关系数一般较高，则可以认为是存在着严重的多重共线性

2、方差膨胀因子

在回归中我们用VIF表示方差膨胀因子

表达式 VIF=1/(1-R^2)

随着多重共线性的严重程度增强，方差膨胀因子会逐渐的变大，一般的当VIF>=10的时候，我们就可以认为存在严重多重共线性；

在R语言中car包中的vif()函数可以帮我们算出这个方差膨胀一找你

这就介绍这两个了，其实还有好多方法，大家可以可以私底下查，或者和我一起交流；

多重共线性的解决办法

因为存在多重共线性，我们还是拟合模型的；当然会有解决办法，这里我就介绍一下常用的方法岭回归；其他的方法也有，这里就不说了；

这里就说说大概的思想，具体推导的步骤这里就不写，有兴趣的可以网上查查；在多重共线性十分严重下，两个共线变量的系数之间的二维联合分布是一个山岭曲面，曲面上的每一个点对应一种残差平方和，点的位置越高，相应的残差平方和越小。因此山岭最高点和残差平方和的最小值相对应，相应的参数值便是参数的最小二乘法估计值，但由于多重共线性的存在最小二乘法估计量已经不适用，一个自然的想法就是应寻找其他的更适合的估计量，这种估计量既要具有较小的方差，又不能使残差平方和过分偏离其极小值。在参数的联合分布曲面上，能满足这种要求的点只能沿着山岭寻找，这就是岭回归法；

这个方法实质是牺牲了无偏性来寻求参数估计的最小方差性；

缺点：通常岭回归方程的R平方值会稍低于普通回归分析，但回归系数的显著性往往明显高于普通回归，在存在共线性问题和病态数据偏多的研究中有较大的实用价值

R语言建模

这里使用可能要使用到car和MASS，由于谢老师已经写了详细的过程，这里我就全程照搬了，偷了个懒，写个代码过程其实也有些累的；

1 分别使用岭回归和Lasso解决薛毅书第279页例6.10的回归问题

cement <- data.frame(X1 = c(7, 1, 11, 11, 7, 11, 3, 1, 2, 21, 1, 11, 10), X2 = c(26, 

29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68), X3 = c(6, 15, 8, 8, 6, 

9, 17, 22, 18, 4, 23, 9, 8), X4 = c(60, 52, 20, 47, 33, 22, 6, 44, 22, 26, 

34, 12, 12), Y = c(78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7, 72.5, 93.1, 

115.9, 83.8, 113.3, 109.4))

cement

##    X1 X2 X3 X4     Y

## 1   7 26  6 60  78.5

## 2   1 29 15 52  74.3

## 3  11 56  8 20 104.3

## 4  11 31  8 47  87.6

## 5   7 52  6 33  95.9

## 6  11 55  9 22 109.2

## 7   3 71 17  6 102.7

## 8   1 31 22 44  72.5

## 9   2 54 18 22  93.1

## 10 21 47  4 26 115.9

## 11  1 40 23 34  83.8

## 12 11 66  9 12 113.3

## 13 10 68  8 12 109.4

lm.sol <- lm(Y ~ ., data = cement)

summary(lm.sol)

## 

## Call:

## lm(formula = Y ~ ., data = cement)

## 

## Residuals:

##    Min     1Q Median     3Q    Max 

## -3.175 -1.671  0.251  1.378  3.925 

## 

## Coefficients:

##             Estimate Std. Error t value Pr(>|t|)  

## (Intercept)   62.405     70.071    0.89    0.399  

## X1             1.551      0.745    2.08    0.071 .

## X2             0.510      0.724    0.70    0.501  

## X3             0.102      0.755    0.14    0.896  

## X4            -0.144      0.709   -0.20    0.844  

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Residual standard error: 2.45 on 8 degrees of freedom

## Multiple R-squared:  0.982,  Adjusted R-squared:  0.974 

## F-statistic:  111 on 4 and 8 DF,  p-value: 4.76e-07

# 从结果看，截距和自变量的相关系数均不显著。

# 利用car包中的vif（）函数查看各自变量间的共线情况

library(car)

vif(lm.sol)

##     X1     X2     X3     X4 

##  38.50 254.42  46.87 282.51

# 从结果看，各自变量的VIF值都超过10，存在多重共线性，其中，X2与X4的VIF值均超过200.

plot(X2 ~ X4, col = "red", data = cement)

接下来，利用MASS包中的函数lm.ridge()来实现岭回归。下面的计算试了151个lambda值，最后选取了使得广义交叉验证GCV最小的那个。

library(MASS)

## 

## Attaching package: 'MASS'

## 

## The following object is masked _by_ '.GlobalEnv':

## 

##     cement

ridge.sol <- lm.ridge(Y ~ ., lambda = seq(0, 150, length = 151), data = cement, 

model = TRUE)

names(ridge.sol)  # 变量名字

## [1] "coef"   "scales" "Inter"  "lambda" "ym"     "xm"     "GCV"    "kHKB"  

## [9] "kLW"

ridge.sol$lambda[which.min(ridge.sol$GCV)]  ##找到GCV最小时的lambdaGCV

## [1] 1

ridge.sol$coef[which.min(ridge.sol$GCV)]  ##找到GCV最小时对应的系数

## [1] 7.627

par(mfrow = c(1, 2))

# 画出图形，并作出lambdaGCV取最小值时的那条竖直线

matplot(ridge.sol$lambda, t(ridge.sol$coef), xlab = expression(lamdba), ylab = "Cofficients", 

type = "l", lty = 1:20)

abline(v = ridge.sol$lambda[which.min(ridge.sol$GCV)])

# 下面的语句绘出lambda同GCV之间关系的图形

plot(ridge.sol$lambda, ridge.sol$GCV, type = "l", xlab = expression(lambda), 

ylab = expression(beta))

abline(v = ridge.sol$lambda[which.min(ridge.sol$GCV)])

par(mfrow = c(1, 1))

# 从上图看，lambda的选择并不是那么重要，只要不离lambda=0太近就没有多大差别。

# 下面利用ridge包中的linearRidge()函数进行自动选择岭回归参数

library(ridge)

mod <- linearRidge(Y ~ ., data = cement)

summary(mod)

## 

## Call:

## linearRidge(formula = Y ~ ., data = cement)

## 

## 

## Coefficients:

##             Estimate Scaled estimate Std. Error (scaled) t value (scaled)

## (Intercept)   83.704              NA                  NA               NA

## X1             1.292          26.332               3.672             7.17

## X2             0.298          16.046               3.988             4.02

## X3            -0.148          -3.279               3.598             0.91

## X4            -0.351         -20.329               3.996             5.09

##             Pr(>|t|)    

## (Intercept)       NA    

## X1           7.5e-13 ***

## X2           5.7e-05 ***

## X3              0.36    

## X4           3.6e-07 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 

## Ridge parameter: 0.01473, chosen automatically, computed using 2 PCs

## 

## Degrees of freedom: model 3.01 , variance 2.84 , residual 3.18

# 从模型运行结果看，测岭回归参数值为0.0147，各自变量的系数显著想明显提高（除了X3仍不显著）

最后，利用Lasso回归解决共线性问题

library(lars)

## Loaded lars 1.2

x = as.matrix(cement[, 1:4])

y = as.matrix(cement[, 5])

(laa = lars(x, y, type = "lar"))  #lars函数值用于矩阵型数据

## 

## Call:

## lars(x = x, y = y, type = "lar")

## R-squared: 0.982 

## Sequence of LAR moves:

##      X4 X1 X2 X3

## Var   4  1  2  3

## Step  1  2  3  4

# 由此可见，LASSO的变量选择依次是X4，X1，X2，X3

plot(laa)  #绘出图数据分析培训

summary(laa)  #给出Cp值

## LARS/LAR

## Call: lars(x = x, y = y, type = "lar")

##   Df  Rss     Cp

## 0  1 2716 442.92

## 1  2 2219 361.95

## 2  3 1918 313.50

## 3  4   48   3.02

## 4  5   48   5.00

# 根据课上对Cp含义的解释（衡量多重共线性，其值越小越好），我们取到第3步，使得Cp值最小，也就是选择X4，X1，X2这三个变量