各类统计方法R语言实现（七）

今天是各类统计方法R语言实现的第七期，我们主要介绍多重共线性、异常观察值的分析和回归模型改进措施。

多重共线性

多重共线性是指线性回归模型中的解释变量之间由于存在强相关关系而使模型估计失真或难以估计准确，它会导致模型参数的置信区间过大，使参数解释较困难。

多重共线性可用VIF（Variance Inflation Factor，方差膨胀因子）进行检测，该指标的经验判断方法：VIF在5到10之间：中度共线性。VIF大于10：重度共线性。

多重共线性解决方法

手动移除出共线性的自变量
逐步回归法
增加样本容量
岭回归或lasso回归
利用因子分析合并变量

#模型拟合
fit<-lm(mpg~hp+wt,data=mtcars)

##展示模型
summary(fit)

## 
## Call:
## lm(formula = mpg ~ hp + wt, data = mtcars)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.941 -1.600 -0.182  1.050  5.854 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 37.22727    1.59879  23.285  < 2e-16 ***
## hp          -0.03177    0.00903  -3.519  0.00145 ** 
## wt          -3.87783    0.63273  -6.129 1.12e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.593 on 29 degrees of freedom
## Multiple R-squared:  0.8268, Adjusted R-squared:  0.8148 
## F-statistic: 69.21 on 2 and 29 DF,  p-value: 9.109e-12

#计算vif
library(car)

## Loading required package: carData

vif(fit)

##       hp       wt 
## 1.766625 1.766625

不存在多重共线性

异常观测值

上次推文已经介绍了异常观测值主要有三类：离群值点、高杠杆值点、强影响点，具体如下：

a.离群点：拟合回归模型对其预测效果不佳（即残差的绝对值较大）。

b.有高杠杆值的变量表明它是一个异常的自变量组合。

c.强影响点表明他对模型参数的估计产生的影响过大。

离群点

之前已经介绍在标准化残差的QQ图中，偏离其他值的异常点可能是离群点，一般认为标准化残差绝对值大于2的点为离群点。

接下来介绍另一种判断离群值的方法，即使用car包中的outlierTest()函数。

fit2<-lm(weight ~ height + I(height^2),data = women)
summary(fit2)

## 
## Call:
## lm(formula = weight ~ height + I(height^2), data = women)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50941 -0.29611 -0.00941  0.28615  0.59706 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 261.87818   25.19677  10.393 2.36e-07 ***
## height       -7.34832    0.77769  -9.449 6.58e-07 ***
## I(height^2)   0.08306    0.00598  13.891 9.32e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3841 on 12 degrees of freedom
## Multiple R-squared:  0.9995, Adjusted R-squared:  0.9994 
## F-statistic: 1.139e+04 on 2 and 12 DF,  p-value: < 2.2e-16

library(car)
outlierTest(fit2)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 15 2.575781           0.025783      0.38675

可以看出15号点p=0.38675，不显著，表明没有离群点。注意：outlierTest（）函数是根据单个最大残差值(绝对值）的显著性来判断是否有离群点，若不显著，则说明数据集中没有离群点，若显著，则必须删除该离群点，然后再检验是否还有其他离群点存在。

高杠杆值点

有高杠杆值的变量表明它是一个异常的自变量组合，即由许多异常的自变量组合起来的异常点，与因变量值没有关系。

高杠杆值的观测点可通过帽子统计量（hat statistic）判断。对于一个给定的数据集，帽子均值为p/n，其中p是模型估计的参数数目（包含截距项），n是样本量。

一般来说，若观测点的帽子值大于帽子均值的2或3倍，则可认定为高杠杆值点。

hat.plot<-function(fit){
  p<-length(coefficients(fit))
  n<-length(fitted(fit))
  plot(hatvalues(fit),main="Index Plot of Hat Values",ylim=c(0,3*p/n)+0.2)

  abline(h=c(2,3)*p/n,col="red",lty=2)
  identify(1:n,hatvalues(fit),names(hatvalues(fit)))
}
hat.plot(fit2)

image.png

## integer(0)

简化代码

##简化代码
hat<-hatvalues(fit2)
hat_mean<-mean(hat)
plot(hat,ylim=c(0,3*hat_mean))
abline(h=c(2,3)*mean(hatvalues(fit2)) , col="red",lty=2)

image.png

水平的两根红线表示帽子均值的2和3倍,可以看出1号点和15号点超过了2倍但没到3倍。

强影响点

表明某点对模型参数的估计产生的影响过大，即移除该点，模型会发生巨大的变化。

检测方法：

Cook距离，或称为D统计量：Cook’s D值大于4/(n-k-1)，则表明它是强影响点，其中n为样本量大小，k是预测变量数目（有助于鉴别强影响点，但并不提供关于这些点如何影响模型的信息）

变量添加图（added variable plot）：对于每个预测变量Xk，绘制Xk在其他k-1个预测变量上回归的残差值相对于响应变量在其他k-1个预测变量上回归的残差值的关系图。

#Cook距离
cutoff<-4/(nrow(women)-length(fit2$coefficients)-2)
plot(fit2,which=4,cook.levels=cutoff)
abline(h=cutoff,lty=2,col="red")

image.png

红线表示4/(n-k-1)，可以发现15号cook距离最大，与上次推文结果一致。

#变量添加图
library(car)
avPlots(fit2,ask=FALSE,onepage=TRUE,id.method="identify")

## Warning in plot.window(...): "onepage" is not a graphical parameter

## Warning in plot.window(...): "id.method" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "onepage" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "id.method" is not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "onepage" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
## a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "onepage" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
## a graphical parameter

## Warning in box(...): "onepage" is not a graphical parameter

## Warning in box(...): "id.method" is not a graphical parameter

## Warning in title(...): "onepage" is not a graphical parameter

## Warning in title(...): "id.method" is not a graphical parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "onepage" is not a
## graphical parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "id.method" is not a
## graphical parameter

## Warning in plot.window(...): "onepage" is not a graphical parameter

## Warning in plot.window(...): "id.method" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "onepage" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "id.method" is not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "onepage" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
## a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "onepage" is not a
## graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
## a graphical parameter

## Warning in box(...): "onepage" is not a graphical parameter

## Warning in box(...): "id.method" is not a graphical parameter

## Warning in title(...): "onepage" is not a graphical parameter

## Warning in title(...): "id.method" is not a graphical parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "onepage" is not a
## graphical parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "id.method" is not a
## graphical parameter

image.png

对于此图，可以想象去掉某一个点之后，直线拟合是否会有大范围变动，此处15号点的影响在所有点中算是比较大的了。

结果整合

car包中的influencePlot（）函数

hat<-hatvalues(fit2)
hat_mean<-mean(hat)

library(car)
influencePlot(fit2,id.method="identify",main="Influence Plot",
              sub="Circle size if proportional to Cook's distance",
              xlim=c(0,3*hat_mean))

## Warning in plot.window(...): "id.method" is not a graphical parameter

## Warning in plot.xy(xy, type, ...): "id.method" is not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
## a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "id.method" is not
## a graphical parameter

## Warning in box(...): "id.method" is not a graphical parameter

## Warning in title(...): "id.method" is not a graphical parameter

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "id.method" is not a
## graphical parameter

image.png

##       StudRes       Hat      CookD
## 1  -0.3527249 0.4647059 0.03883656
## 2  -1.5156988 0.2680672 0.25310078
## 13 -1.5312900 0.1656755 0.13956756
## 15  2.5757809 0.4647059 1.30646190

本质上是将三个值绘制在一张图里。

纵坐标超过2或小于-2的点可被认为是离群点，水平轴超过2倍或3倍帽子值均值的点有高杠杆值。圆圈大小与影响成比例，圆圈很大的点可能是对模型估计造成的不成比例影响的强影响点。

回归模型改进措施

主要有四种方法：

（1）删除异常值.

（2）变量变换。

（3）添加或删除变量。

（4）使用其他回归方法。

删除异常值

通常删除离群点和强影响点，直到拟合较满意。

当然删除要谨慎，可以探究产生异常值的原因。

变量变换

可以尝试各类变换方法，使变量满足正态性、线性、同方差性假设，可以尝试之前各类统计方法R语言实现（四）介绍的方法，但是变量变换之后需要有具体意义。

添加或删除变量

可以删除多重共线性的变量（根据VIF），也可以岭回归或lasso回归。

使用其他回归方法

存在离群点或强影响点，可使用稳健回归代替最小二乘回归。

违背了正态性假设，可以用非参数回归模型。

存在显著非线性，可以使用非线性模型。

违背了误差独立性假设，可以使用专门研究误差结构的模型，如时间序列模型或多层次回归模型。

最后，还能依据数据的分布形式选择不同的广义线性模型。

好了，今天的R语言实现统计方法系列推文暂时告一段落，我们下次再见吧！小伙伴们如果有什么统计上的问题，或者如果想要学习什么方面的生物信息内容，可以在微信群或者知识星球提问，没准哪天的推文就是专门解答你的问题哦！