R语言_OLS回归分析

R语言_学习笔记

Content

(1). Part 1. 回归简介-------基础回顾、OLS适用环境
(2). Part 2. 具体方法-------简单线性回归、多项式回归、多元线性回归
(3). Part 3. 回归诊断-------线性模型的验证、多重共线性、异常观测值、改进方法
(4). Part 4. 改进措施-------删除观测点、变量变换、增删变量、尝试其它方法
(5). Part 5. 选择“最佳”的回归模型---模型比较、变量选择
(6). Part 6. 深层次分析----交叉验证、预测变量的相对重要性

Part 1. 回归简介

1、基础回顾

• 定义：回归分析是指用一个或多个解释变量来预测响应变量。
• OLS（普通最小二乘）：包括简单线性回归、多项式回归、多远线性回归
• 漂亮的OLS模型应满足：
  
  • 正态性 ： 对于固定的自变量值，因变量值成正态分布。
  • 独立性 ： 因变量值（残差）之间相互独立。
  • 线性 ： 因变量与自变量之间为线性相关。
  • 同方差性 ：因变量的方差不随自变量的水平不同而变化。
    （如何判断，请参见Part.3）

2、OLS适用环境

• OLS回归是通过预测变量的加权和来预测量化的因变量，其中权重是通过数据估计而得的参数。
• OLS擅于处理的问题类型：

 铺路表面的面积与表面盐度有什么关系？
 哪些因素可以解释各地的啤酒价格差异？
 血压、盐摄入量与年龄的关系是怎样的？对于男性和女性是相同的吗？
 ......

Part 2. 具体方法

1、简单线性回归

> library(MASS)
> attach(Cars93)
> par(bty="l")    #坐标轴呈现“L”状
> fit1 <- lm(RPM~Horsepower)  
> plot(Horsepower,RPM,
       xlab="Horsepower",
       ylab="RPM" )
> abline(fit1,col="darkorange2",lwd=2)
> summary(fit1)
Call:
lm(formula = RPM ~ Horsepower)
Residuals:
    Min      1Q  Median      3Q     Max 
-1479.0  -462.3   -54.8   487.7  1226.4 
Coefficients:
            Estimate Std. Error t value   Pr(>|t|)    
(Intercept) 5220.524    182.581   28.59   <2e-16 ***
Horsepower     0.418      1.194    0.35    0.727    
Residual standard error: 599.6 on 91 degrees of freedom
Multiple R-squared:  0.001346,  Adjusted R-squared:  -0.009628 
F-statistic: 0.1227 on 1 and 91 DF,  p-value: 0.727

2、多项式回归

• 从上面的图，我们可以看到模型的拟合效果不好。现用高大上的scatterplot()函数来直观说明这一问题。

> library(car)
> library(MASS)
> attach(Cars93)
> scatterplot(RPM~Horsepower,
              spread=FALSE,
              lty=2,
              pch=19,
              lwd=1.6,
              xlab="Horsepower",
              ylab="RPM",
              #smoother=gamLine,
              span=0.5,       
              #span越大，产生的函数就越平滑
              col=c("blue", "darkorange2","olivedrab3") )
  #此处默认的平滑方法是LOESS（局部拟合方法），可以设置参数
  #smoother=gamLine 调用另一种平滑方法。

实线为平滑拟合，虚线为线性拟合。

• R表达式中常用的符号：

：   表示预测变量的交互项
*    表示所有可能交互项的简洁表示
     y~x*w*z 可展开为 y~x+w+z+x:w+x:z+w:z+x:w:z
^    表示交互项达到某个次数  
     y~(x+y+z)^3 可展开为 y~x+w+z+x:w+x:z+w:z
.    表示包含除因变量的所有变量的简洁表示
     y~.表示y~x+w+z  （假设数据框中只有有x,y,z,w四个变量） 
+0   删除截距项
I()  从算术的角度来解释括号中的元素
     y~x+I((z+w)^2) 可展开为 y~x+h,h是z和w的平方和

• 现利用多项式回归重构模型

> library(car)
> library(MASS)
> attach(Cars93)
> fit2 <- lm(RPM~Horsepower + I(Horsepower^3))
> plot(Horsepower,RPM,
       main="RPM&Horsepower",
       xlab="Horsepower",
       ylab="RPM",
       col="olivedrab3",
       lwd=2
      ) 
> lines(sort(Horsepower),sort(fitted(fit2)),type="l",col="darkorange2",lwd=2)
由于是电脑作图，其思路和人类有所不同。所以一定要现对自变量和因变量进行排序，否则画出的图会很奇怪。

添加多项式后，拟合程度有所改善，但依旧不够好。原因在于预测变量太少，无法准确衡量响应变量。

3、多元线性回归

• 做多元回归分析的第一步，是先对变量做相关性检验(cor.test)。

> library(MASS)
> library(psych)
> attach(Cars93)
> car <- Cars93[,c("Weight","Horsepower","RPM","EngineSize")]
> corr.test(car,method="kendall")
#method参数还可以选择默认的“pearson”或者“spearman”
Call:corr.test(x = car, method = "kendall")
Correlation matrix     #相关系数矩阵
           Weight     Horsepower     RPM     EngineSize
Weight       1.00       0.63        -0.31       0.74
Horsepower   0.63       1.00        -0.05       0.65
RPM         -0.31      -0.05         1.00      -0.40
EngineSize   0.74       0.65        -0.40       1.00
Sample Size 
[1] 93
Probability values (Entries above the diagonal are adjusted for multiple tests.) 
           Weight     Horsepower     RPM     EngineSize
Weight          0       0.00         0.01       0
Horsepower      0       0.00         0.63       0
RPM             0       0.63         0.00       0
EngineSize      0       0.00         0.00       0
从结果看出，RPM与Horsepower的相关系数为-0.05，p值为0.63,并非显著地不为0。

• 多元线性回归

> library(MASS)
> attach(Cars93)
> fit3 <- lm(RPM~Horsepower+Weight+EngineSize)
> summary(fit3)
Call:
lm(formula = RPM ~ Horsepower + Weight)
Residuals:
     Min       1Q   Median       3Q      Max 
-1351.75  -284.51    29.87   320.96   879.82 
Coefficients:
             Estimate Std. Error  t value   Pr(>|t|)    
(Intercept) 7121.9417   257.1840   27.692   < 2e-16 ***
Horsepower     8.8515     1.3129    6.742   1.45e-09 ***
Weight        -1.0135     0.1166   -8.695   1.49e-13 ***
Residual standard error: 444.5 on 90 degrees of freedom
Multiple R-squared:  0.4572,    Adjusted R-squared:  0.4452 
F-statistic: 37.91 on 2 and 90 DF,  p-value: 1.141e-12

虽然各系数的都显著不为0,但是R-squared的值不够优秀。
这可能是由于没有考虑到交互项的影响导致的。

• 考虑了交互项的多元线性回归

> library(MASS)
> attach(Cars93)
> fit4 <- lm(RPM~Horsepower+Weight+EngineSize+
             Horsepower:Weight:EngineSize)
> summary(fit4)
Call:
lm(formula = RPM ~ Horsepower + Weight + EngineSize + Horsepower:Weight:EngineSize)
Residuals:
    Min      1Q  Median      3Q     Max 
-825.51 -172.25   42.86  209.36  758.08 
Coefficients:
                               Estimate Std. Error t value Pr(>|t|)    
(Intercept)                   6.494e+03  2.652e+02  24.482  < 2e-16 ***
Horsepower                    9.088e+00  1.595e+00   5.697 1.59e-07 ***
Weight                       -2.404e-01  1.127e-01  -2.132   0.0358 *  
EngineSize                   -8.001e+02  1.182e+02  -6.769 1.39e-09 ***
Horsepower:Weight:EngineSize  2.459e-04  1.229e-04   2.000   0.0486 *  
Residual standard error: 305.8 on 88 degrees of freedom
Multiple R-squared:  0.7487,    Adjusted R-squared:  0.7373 
F-statistic: 65.56 on 4 and 88 DF,  p-value: < 2.2e-16

• 添加了交互项后，各项系数仍然显著不为0，模型的R-squared提高到了0.7373。
• 仔细一点可以发现，交互项的p值=0.0486，表明显著性不高，可以考虑删除交互项。
• 暂时把回归模型确定为：

fitc <- lm(RPM~Horsepower+Weight+EngineSize)

Part 3. 回归诊断

1、假设条件诊断
对假设条件再啰嗦一下:
- 正态性 ： 对于固定的自变量值，因变量值成正态分布。
- 独立性 ： 因变量值（残差）之间相互独立。
- 线性 ： 因变量与自变量之间为线性相关。
- 同方差性 ：因变量的方差不随自变量的水平不同而变化。

a、基本诊断方法

对lm()函数返回的对象使用plot()函数，生成用来评价回归模型拟合情况的四幅图形。

> par(mfrow=c(2,2))
> plot(lm(weight~height+I(height^2),data=women))

这四幅图就显示回归模型“weight=height+height^2”非常漂亮。可以从OLS回归的统计假设来理解：

• 正态性
  若满足此假设，残差值也应是一个均值为0的正态分布，正态Q-Q图（Normal Q-Q右上）的点应该落在呈45°角的直线上。

• 独立性
  此假设无法从四幅图中验证。

• 线性
  若满足此假设，则残差图与拟合图（Residuals vs Fitted,左上）应成随机分布。

• 同方差性
  若满足此假设，则位置尺度图（Scale-Location Graph，左下）中，水平线周围的点应随机分布。

• 最后一幅“残差与杠杆图”（Residuals vs Leverge,右下）,可用来鉴别出“离群点”、“高杠杆值点”、“强影响点”

  • 离群点。如果模型中有一个观测点是离群点，则表明该模型对其预测效果不佳（产生了巨大的或正或负的残差）
  • 高杠杆值点。如果一个预测值有很高的杠杆值，则表明它是一个离群点。
  • 强影响点。表明它对模型参数的估计产生的影响过大。强影响点可以通过Cook距离，即Cook'sD统计量来鉴别。

> 如何删除上述不好的点？假设要删除的点是点5,点17(从图中看出)
> newfit <- lm(weight~height+I(height^2),data=women[-c(5,17),])

上述的诊断方法太弱，使用起来也很模糊。下面介绍来自car包的更高大上的诊断方法

2、高大上的诊断方法

• 正态性

> library(car)
> library(MASS)
> attach(Cars93)
> fitc <- lm(RPM~Horsepower+Weight+EngineSize)
> qqPlot(fitc,
         labels=row.names(Cars93)，  #方便在交互模式下点击的点，可以显示出其name或number。
         id.method="identify",
         simulate=TRUE,main="Q-Q Polt")
#id.method="identify",设置当图绘制完成后，用鼠标单击图形内的点，将会标注函数中label选项中的设定值

除了57号点之外，所有的点都离直线很近，并都落在置信区间内，表明正态性假设符合得很好。

现在重点考察57号点：

> car["57",]
    Weight    Horsepower      RPM      EngineSize
57    2895        255         6500          1.3
> fitted(fitc)["57"]
      57 
7447.945 
可以看到，57号点的RPM为6500,而模型预测的RPM为7447.945。
这可能就要考虑其它的影响因素。

• 误差独立性

> durbinWatsonTest(fitc)
 lag Autocorrelation D-W Statistic p-value
   1      -0.1076394      2.160802   0.398
 Alternative hypothesis: rho != 0

p=0.398 >0.05，说明无自相关性，误差项之间独立。
lag=1，表明数据集中，每个数据都是与其后一个数据进行比较。

• 线性

> crPlots(fitc)

crPlot()函数绘制的成分残差图表明，线性模型形式对该数据集看似是合适的。

• 同方差性

> ncvTest(fitc)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 15.54264    Df = 1     p = 8.066492e-05 

当p>0.05时，表明满足同方差性假设。此时的模型还不满足该假设。

• 综合验证

> library(gvlma)
> gvmodel <- gvlma(fitc)
> summary(gvmodel)
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 
Call:
 gvlma(x = fitc) 
                     Value     p-value        Decision
Global Stat         11.0697    0.02579    Assumptions NOT satisfied!
Skewness            2.9119     0.08793    Assumptions acceptable.
Kurtosis            1.2877     0.25648    Assumptions acceptable.
Link Function       6.7566     0.00934    Assumptions NOT satisfied!
Heteroscedasticity  0.1134     0.73627    Assumptions acceptable.

从输出结果看到，“Global Stat”的p值小于0.05，所以不通过综合检验。

3、多重共线性
是指回归模型中，由于自变量之间存在较强的相关关系，而使模型估计失真。

vif(fitc)
Horsepower     Weight     EngineSize 
  2.418150    3.925249     3.842155 

vif的值小于4,表明不存在多重共线性。
如果vif的值都大于4,应使用岭回归。

4、异常观测值

• 离群点

> outlierTest(fitc)
       rstudent     unadjusted p-value    Bonferonni p
57     -3.781735         0.00028332         0.026349

从输出结果来看，utlierTest()总是找出残差绝对值最大的点。只有当Bonferonni p<0.05时，才为离群点。例如此例中的57号点被认定为离群点。

• 高杠杆值点

> hat.plot <- function(fit){
    p <- length(coefficients(fit))
    n <- length(fitted(fit))
    plot(hatvalues(fit),main="Index Plot of Hat Values")
    abline(h=c(2,3)*p/n,col="red",lty=2)
    identify(1:n,hatvalues(fit),names(hatvalues(fit)))
 }
> hat.plot(fitc)
#高于两条虚线的点，都可以被认为是高杠杆值点。

• 强影响点

使用Cook距离（D统计量）寻找强影响点,优点是可以直接判别。

> cutoff <- 4/(nrow(car)-length(fitc$coefficients)-2)  
#这是判断条件，其中nrow(car)表示样本容量，length(fitc$coefficients)表示预测变量数目
> plot(fitc,which=4,cook.levels=cutoff)
> abline(h=cutoff,lty=2,lwd=2,col="orange")

• 将三种异常观测值绘于一张图上

> library(car)
> library(MASS)
> attach(Cars93)
> fitc <- lm(RPM~Horsepower+Weight+EngineSize)
> influencePlot(fitc,id.method="identify",main="Influence Plot")
> abline(h=-2,col="light blue",lwd=4)
> abline(h=2,col="light blue",lwd=4)
> abline(v=0.3,col="orange",lwd=4)

• 在浅蓝色线之外的点，是离群点;
• 在橙色线右边的点，是高杠杆值点;
• 圆圈很大的点，可能是对模型参数的估计造成不成比例影响的强影响点;

Part 4. 改进措施

• 删除观测点
  针对之前的回归模型，现删除离群点和强影响点后，重新评测模型的拟合度。

> library(car)
> library(MASS)
> library(gvlma)
> fitc2 <- lm(RPM~Horsepower+Weight+EngineSize,
              data=Cars93[-c(11,48,57,19,89),])
> qqPlot(fitc2,labels=row.names(Cars93),
         id.method="identify",
         simulate=TRUE,
         main="Q-Q Polt" )

• 变量变换

使用boxTidwell()函数可以寻找到合适的自变量的次幂，来改善模型的线性关系。

> library(car)
> library(MASS)
> attach(Cars93)
> boxTidwell(RPM~Horsepower+Weight+EngineSize)
           Score Statistic     p-value    MLE of lambda
Horsepower       -2.860639    0.0042279     0.0347044
Weight            1.128053    0.2592974    -0.9517346
EngineSize        1.299718    0.1936977     0.5270448

若P值<0.05,则表明有必要添加对应的幂到模型中。
此处应将 “0.0347044” 添加为Horsepower的指数。
综上，模型的最终版确定为：

fitchloe <- lm(RPM~I(Horsepower^0.0347044)+Weight+EngineSize,
              data=Cars93[-c(11,48,57,19,89),])

对最终模型进行评测：

> library(car)
> library(MASS)
> library(gvlma)
> fitchloe <- lm(RPM~I(Horsepower^0.0347044)+Weight+EngineSize,
                 data=Cars93[-c(11,48,57,19,89),])
> qqPlot(fitchloe,labels=row.names(Cars93),
         id.method="identify",
         simulate=TRUE,
         main="Q-Q Polt" )

通过正态性假设检验

> durbinWatsonTest(fitchloe)
 lag Autocorrelation D-W Statistic p-value
   1     -0.06875225      2.105465    0.63
 Alternative hypothesis: rho != 0

通过残差独立性假设检验

> crPlots(fitchloe)

通过线性假设检验

> ncvTest(fitchloe)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 2.558461    Df = 1     p = 0.1097053 

通过同方差性检验

> gvmodel <- gvlma(fitchloe)
> summary(gvmodel)
Call:
lm(formula = RPM ~ I(Horsepower^0.0347044) + Weight + EngineSize, 
    data = Cars93[-c(11, 48, 57, 19, 89),])
Residuals:
    Min      1Q  Median      3Q     Max 
-852.46 -186.74  -23.39  166.94  845.95 
Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)             -4.656e+04  4.406e+03 -10.567  < 2e-16 ***
I(Horsepower^0.0347044)  4.639e+04  3.892e+03  11.920  < 2e-16 ***
Weight                  -5.511e-01  1.183e-01  -4.661  1.1e-05 ***
EngineSize              -5.494e+02  6.082e+01  -9.034  3.2e-14 ***
---
Residual standard error: 314 on 89 degrees of freedom
Multiple R-squared:  0.7321,    Adjusted R-squared:  0.7231 
F-statistic: 81.08 on 3 and 89 DF,  p-value: < 2.2e-16
Call:
 gvlma(x = fitchloe) 
                       Value    p-value                Decision
Global Stat        4.3792938    0.35710        Assumptions acceptable.
Skewness           0.0004179    0.98369        Assumptions acceptable.
Kurtosis           1.3559167    0.24425        Assumptions acceptable.
Link Function      2.9540798    0.08566        Assumptions acceptable.
Heteroscedasticity 0.0688794    0.79298        Assumptions acceptable.

通过所有检验，各参数均符合显著性条件

> vif(fitchloe)
I(Horsepower^0.0347044)       Weight                EngineSize 
     3.058320                3.540496                3.713773 

且模型不存在多重共线性

• 增删变量
  主要用于处理多重共线性，方法有两种：
  
  • 删除√vif>2的变量;
  • 岭回归，此法高大上，值得关注;

Part 5. 选择“最佳”的回归模型

1、模型比较

> fitchloe <- lm(RPM~I(Horsepower^0.0347044)+Weight+EngineSize,data=Cars93[-c(11,48,57,19,89),])
> fitcx <- lm(RPM~I(Horsepower^0.0347044)+EngineSize,data=Cars93[-c(11,48,57,19,89),])
> anova(fitchloe,fitcx)
Analysis of Variance Table
Model 1: RPM ~ I(Horsepower^0.0347044) + Weight + EngineSize
Model 2: RPM ~ I(Horsepower^0.0347044) + EngineSize
  Res.Df      RSS Df Sum of Sq      F    Pr(>F)    
1     89  8775711                                  
2     90 10917425 -1  -2141714 21.721 1.101e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#由于p<0.05，表明应该将Weight作为自变量添加到模型中。

> AIC(fitchloe,fitcx)
          df      AIC
fitchloe  5    1339.228
fitcx     4    1357.537

fitchloe模型的AIC值更小，所以应该优先选择该模型。

2、变量选择
在原始数据集中，如何挑选用于解释因变量的自变量？？

• 较弱的逐步回归

> library(MASS)
> library(leaps)
> library(car)
> fitc44 <- lm(RPM~Horsepower+Weight+EngineSize+Price+MPG.city+MPG.highway+Rev.per.mile
              +Fuel.tank.capacity+Length+Width+Turn.circle+Rear.seat.room)
> stepAIC(fitc44,direction="backward")

最后的结果：

Step:  AIC=1029.84
RPM ~ Horsepower + Weight + EngineSize + Length + Rear.seat.room
                 Df     Sum of Sq        RSS       AIC
<none>                                 6556343    1029.8
- Length          1      151055        6707398    1029.9
- Rear.seat.room  1      482056        7038399    1034.3
- Weight          1      1397774       7954117    1045.4
- EngineSize      1      8200814       14757157   1101.7
- Horsepower      1      14576782      21133125   1134.3
Call:
lm(formula = RPM ~ Horsepower + Weight + EngineSize + Length + 
                   Rear.seat.room)     #逐步回归算法的推荐模型。
Coefficients:
(Intercept)  Horsepower  Weight  EngineSize  Length  Rear.seat.room 
 4840.5086   14.8568    -0.5188  -730.1240   5.3970    31.4344  

• 高大上的全子集回归

> fitc4 <- regsubsets(RPM~Horsepower+Weight+EngineSize+Price+MPG.city+MPG.highway+Rev.per.mile
                          +Fuel.tank.capacity+Length+Width+Turn.circle+Rear.seat.room)
> subsets(fitc4,statistic="cp",main="Cp Plot for ALL Subsets Regression")
#subsets()函数的作用是找到最佳的子集预测。
#statistic有五个不同的值，"bic":贝叶斯信息标准; "cp":Mallow的Cp;"adjr2":R^2调整自由度;"rsq":未经调整的R^2;"rss"，残差平方和; 更多时候，选择cp的精确度较高，且无偏倚。
> abline(1,1,lty=2,col="red")

注意在保存图形时，要保存为eps格式的大尺寸矢量图。
此处受到环境限制，只能上传小尺寸的图。

Part 6. 深层次分析

1、交叉验证

• “交叉验证”：将一定比例的数据挑选出来作为训练样本，另外的样本作保留样本，现在训练样本上获取回归方程，然后在保留样本上做预测。
• 利用shrinkage()函数做R平方的k重交叉验证。

> shrinkage(fitchloe)
Orginal R-square =  0.7321225 
10 Fold Cross-Validated R-square = 0.7084592 
Change = 0.02366325 

分析结果显示，基于初始样本的R^2等于0.7321225，很漂亮。交叉验证后的R^2等于0.7084592，变化不大。这说明回归模型fitchloe有很好的预测能力。
（由于观测值是被随机分配到k个群组中，所以每次运行shrinkage()函数，得到的结果会有少许不同。）

2、预测变量的相对重要性
用高大上的relweights()函数分析回归模型中各个预测变量的相对权重。

> relweights(fitchloe,col="cadetblue1")
                            Weights
I(Horsepower^0.0347044)     27.45912
Weight                      28.13227
EngineSize                  44.40862