R语言实战第八章

这一张主要介绍了回归分析，

  1. 回归的类型：回归包括哪些回归，比如一元线性回归、多项式回归、多元线性回归、Logist回归、泊松回归、多变量回归等，这张主要以一元、多项式和多元线性回归举例
  2. 因为你做回归的模型都是有统计假设的，而你做了这个模型，就要确保这个回归方程满足回归模型的统计假设。
  3. 如果有不满足回归假设的情况，可以通过添加变量、变量转换、删除离群点、强影响点来使其满足回归模型的统计假设。
  4. 得到回归模型后，如果判断回归模型的好坏，比如多个回归模型，如何比较哪个更好，以及选择哪些变量
  5. 得到回归模型后，还要评估模型的泛化能力和变量的相对重要性。
    本章所有的例子都是以线性模型举例

1. 回归的类型和示例

image.png

本章的示例是以OLS法通过一系列的预测变量来预测相应变量


image.png

n是观测的数目，k为预测变量的数，i是观测值。
OLS的统计假设要求是：
● 正态性，因变量是
● 独立性，因变量值之间
● 线性，自变量和因变量
● 同方差性，因变量的方差不随自变量的水平不同而变化。称为不变方差。
线性模型通过lm拟合，fit<-lim(formula,data)
formula的形式有多重多样，最简单的形式是y~X1+X2...


image.png

此外我们得到回归模型后有一系列的函数可以得到关于模型的信息
image.png
  1. 一元线性模型
    但实际上一开始你是没有办法判定到底该选择什么样子的模型的，你只能去不停地尝试，不停地修正，然后看看哪些最符合这些模型，比如本章一开始的例子是用身高预测体重。
fit <- lm(weight ~ height, data=women)
summary(fit)
plot(women$height,women$weight,xlab="Height (in inches)", ylab="weight (in pounds)")
abline(fit)
image.png

image.png

我们可以看到这个简单的线性模型应该是wight = -87.51667+3.45x height，这个线性模型呢呢狗狗解释0.991的模型情况，F检验也是显著的，而两个参数截距intercept和height也都是显著的。
这里有一个Resudual standard error残差标准差1.525，可认为是模型用身高预测体重的平均误差。F统计量检验的所有的预测变量预测的相应变量是否都在某个几率水平之上。
但是，我们可以看到实际上，通过一条直线并不能很好的拟合所有的数据，这个图形暗示我们可以通过使用一个含有弯曲曲线来提高预测精度。比如通过模型Y=β+β1X+β2X2

  1. 多项式线性回归
fit2<-lm(weight~height+I(height^2),data=women)
summary(fit2)
plot(women$height,women$weight,xlab="Height (in inches)",ylab="Weight (in lbs)")
lines(women$height,fitted(fit2))
image.png

所以新的预测模型就是weight= 261.878-7.3483 x height + 0.08306(height)^2
可以看到二项式的R方更大，其能够更好的解释数据。
针对于二元数据，这里介绍了car包的scatterplot，让更方便的绘制二元关系图


image.png
  1. 多元线性回归。
    这里以一个州的犯罪率与其他因素的关系。
    实际上在进行多元线性回归的时候，我们要看一些各个变量之间的相关系数，然后通过散点图矩阵查看一些他们之间的关系
states<-as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
cor(states)
library(car)
scatterplotMatrix(states,spread=FALSE,main="Scatter Plot Matrix")
image.png

可以看到只在非对角线的位置绘制变量间的散点图，天剑了平滑和线性拟合曲线，在对角线区域是绘制了密度图和轴须图。
有交互项的多元线性回归

fit <- lm(mpg~hp+wt+hp:wt,data=mtcars)
summary(fit)
image.png

我们可以看到hp:wt之间，即马力和车重是显著的，这代表响应变量与其中一个预测变量的关系受到另一个预测变量的影响。也就是说马力与mpg之间的关系实际上还受到车重的影响，车重不同，马力和mpg之间的关系不同。
这里它通过举例不同的车中，2.2,3.2,4.2的数值来看，mpg与hp之间的关系。

fit <- lm(mpg~hp+wt+hp:wt,data=mtcars)
summary(fit)
library(effects)
#plot(effect("hp:wt",fit,list(wt=c(2.2,3.2,4.2))),multiline=TRUE)
plot(effect("hp:wt",fit,,list(wt=c(2.2,3.2,4.2))),multiline=TRUE)
image.png

从这个我们可以看到hp和mpg本身是负相关的，但是随着wt越来越大，这种负相关会逐渐减少，甚至最后当wt=4.2的时候，hp和mpg之间的关系接近与0。
总结：从上面，我们讲解了回归有哪些类型：一元线性、多项式、多元、泊松、Logists等，以一元线性模型举例，公式应该是什么样的，公式有什么样子的形式，包括模型一些常用的函数，比如summary，fitted等函数，而后讲述了一元线性、多项式和多元，在做多元的时候，我们首先需要关注下变量之间的相关性以及可以通过car包的scatterplotMatrix来展示他们，包括多元线性的交互项，以及交互项的展示方式通过effects包的effect函数。

2. 回归模型的统计假设诊断

我们虽然通过模型进行了数据的拟合，但是实际上诸多模型都有很多前置条件，比如上面的OLS回归模型，我们只是单纯的通过模型的参数推断，以及数据的显著性和参数对数据的解释程度来推断这个模型可能符合数据，但是实际还要确定这批数据多大程度上满足数据的统计假设信息。
这里面提供了两种方式检验模型的统计假设诊断，一种是标准方法，第二种是car包的诸多函数。
当然了本章的介绍也是关于线性模型的统计假设诊断。
我在这里重述，线性模型的要求是什么，首先要求因变量满足正态性、方差同质、独立性以及自变量和因变量之间是线性关系。所以关于线性模型的统计诊断也会从这几个方面着手。
也就是说如果你用一元线性回归推断得出的线性方程，其参数能够极大程度解释数据，而且对数据的拟合很好，参数也具有显著性；然后已用的这个线性模型的统计推断也符合，那么这个模型就不错，可以用，但是是不是最好的模型，我们无从得知，参数的数量多少我们也无从得知，模型的泛化能力（预测新数据的能力怎么样我们也无从得知），当然这些我们后面会解释。

2.1 标准方法

fit<-lm(weight~height,data=women)
summary(fit)
par(mfrow=c(2,2))
plot(fit)
image.png

● 正态性，右上角的Q-Q图，因变量成正态性，其残差只应该是均值为0的正态分布，这是标准残差的概率图，如果满足正态性，数据的点应该落在45°直线上。若不是违反正态性。
● 独立性，数据的独立性需要自我进行判断，比如我们的因变量是女性的体重，我们没有理由相信一个女生的体重会影响另一个女生的体重，但是也有例外，比如一个家庭，然后每天吃的一样，但是两个女生，一个吃的太多，一个吃的太少，那么这两个女生的体重肯定不是独立的。
● 线性，坐上的图，如果因变量与自变量线性相关，残差只与拟合值之间应该是没有任何关系，但是这个数据中呈现了一个曲线关系，这表明了这个数据可能不是一个线性相关的。
● 同方差性，左下，满足方差，标准差和拟合值之间应该是没有任何关系的，这个图似乎满足这个条件。
在最后一幅图中，关于杠杆和标准差，实际这个图的目的是鉴别数据的离群点、杠杆点和强影响点。。
离群点表示拟合的不加，代表有很大的残差，强影响点代表做了这个点会极大程度上影响模型的参数；而杠杆点代表的是预测变量的组合。
但实际右下这个图的实用性并不好。
这里我想介绍一下，残差(residuals)、标准差（standardized residuals）和学生化残差（studentized Residuals），其中残差就是数据值和拟合值之间的差异，标准差是将残差标准化，学生化残差更好的找到离群点、强影响点和杠杆值。
然后，重新用二项式进行了拟合查看数据的好坏情况

newfit <- lm(weight~height+I(height^2),data=women)
image.png

我们从数据中可以看到由上图13点不符合正态性，右下角15号店的Cook distance更大，是强影响点。
这可能代表，我们可以将13和15号提出来满足更好的数据拟合，但是删除要慎重，因为你的数据是观测得到的。

2.2 car包一系列函数评估

image.png

image.png

● qqPlot()可以用来看数据的正态性分布
● durbinWatsonTest()查看数据的独立性与否
● crPlots()可以查看自变量与因变量之间是否线性相关
● ncvTest()和spreadLevelPlot()来查看方差的同质性
● outlierTest()用来查看数据的离群点
● avPlotss()通过变量添加图来查看强影响点。
● influencePlot()来直接判断数据的离群点、杠杆点和强影响点。
下面就具体举例说明

  1. 验证数据的正态性
    方法1：qqPlot()进行验证
library(car)
states<-as.data.frame(state.x77[,c("Murder","Population",
                                   "Illiteracy","Income","Frost")])
fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)
qqPlot(fit)
image.png

方法二：通过residplot()绘制血生化残差柱状图，直方图，添加正态性、密度行曲线而后周虚线
实际上，上面普通方法的plot和qqPlot()的假定都是说因变量是符合正态分布的，那么理论上其残差也应该是符合正态分布的，这个时候就可以直接将studentized Residual进行可视化与正态分布对比，如果符合，那么也就是符合正态分布

residplot<-function(fit,nbreaks=10){
  z<-rstudent(fit)
  hist(z,breaks=nbreaks,freq=FALSE,xlab="Studentized Residual",main="Distribution of Errors")
  rug(jitter(z),col="brown")
  curve(dnorm(x,mean=mean(z),sd=sd(z)),add=TRUE,col="blue",lwd=2)
  lines(density(z)$x,density(z)$y,col="red",lwd=2,lty=2)
  legend("topright",legend=c("Normal Curve",""))
}
image.png
  1. 误差独立性，利用durbinWatsonTest()
    因变量是否独立最好还是自身的背景知识，正如我上面提到的女性身高和体重的预测，但实际上R语言也通过car包提供了一个可做Durbin-Watson检验的函数，能够检测误差的序列相关性。
    durbinWatsonTest(fit)
  2. 线性关系
    自变量和因变量之间要求是线性关系，可以通过成分残差图crPlots()函数查看
library(car)
crPlots(fit)
image.png

如果图形存在非线性，就说明预测变量的函数的建模不够充分，需要添加一些曲线成分，比如多项式等等。
通过成分残差图可以看到自变量和因变量之间应该是线性的

  1. 同方差性，car包通过ncvTest()和spreadLevelPlot()两个函数说明


    image.png

    如果点在拟合曲线钟斌随机分布，则说明方差不变，如果看到一个非水平的情况，说明需要进行幂次变换。
    综上，上面通过了两种方式，即基础方法和car包的一系列方法来检测线性模型的正态性、方差独立性、线性关系以及同方差性的情况。
    线性模型统计假设的综合检验
    但是为了更好的进行判断线性模型的统计假设，R也提供一些R包能够一次性检测上面的统计假设，这里用的包是gvlma的gvlma()函数。这个函数对线性模型进行综合检验，此外还能做偏斜度和异方差性的检验。

library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)
image.png

从上面我们可以看到数据基本上满足OLS回归模型的所有统计假设
多重共线性的解释
这个多重共线性我不太理解。


image.png

发现回归整体显著(F检验p < 0.001)
但是DOB和Age的回归系数都不显著
原因：DOB和Age之间高度相关（几乎等价）
所谓的多重共线性，就是两个或多个自变量高度先关（即它们之间重复提供信息），模型就出现了多重共线性。
DOB和Age本质上是反向线性相关的
Age=CurrentYear-YearOfBirth
但是当把DOB和Age同时作为自变量放进模型时：
● 模型整体拟合（F检验）显著说明这两个变量合起来确实对握力有解释力
● 但是个别变量的系数不显著，。
实际上模型在问这样的问题，年龄固定，DOB对握力有没有影响。但实际上这个逻辑不成立(年龄固定，DOB不能变。)所以模型在数学上很难区分是谁真正贡献了解释力，导致每个系数的标准误很大，置信区间很宽，p值不显著。是这样的原因。

2. 3数据中异常点

而实际上统计模型的满足过程中，经常会有一些异常点，但实际上我们期望回归分析能够覆盖这些异常点，包括离群点、高杠杆值点和强影响点，但实际上有一些影响点会严重影响回归模型，而对这些异常点进行剔除会对回归模型造成极大的影响。
这些异常点，包括了离群点、杠杆点和强影响点。

离群点

离群点是预测模型效果不佳的点，通常会具有比较大的残差值，即模型的预测值和实际的观察值之间的差别是很大的。
我们可以通过Q-Q图来查看离群点，而离群点非常粗糙的判断方法是：标准化残差值大于2或者小于-2的点可能是离群点。
car包中的outlierTest()可以查看数据的离群点


image.png

高杠杆值点

高杠杆值点是与预测变量有关的离群点。这些所谓的高杠杆值点是由异常的预测变量值组合起来的。
高杠杆值点可以通过帽子统计量(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")
  abline(h=c(2,3)*p/n,col="red",lty=2)
  #identity(1:n,hatvalues(fit),names(hatvalues(fit)))
}
hat.plot(fit)
image.png

强影响点

强影响点事值哪些对模型参数估计值影响很大的点，比如如果移除掉模型的一个观测点的时候，模型会发生巨大的改变，这个时候需要检测模型中是否有强影响点的。
两种方法检测强影响点：Cook距离和变量添加图。

  1. Cook距离也称为D统计量
    Cook的D值大于4/(n-k-1)表示其是强影响点，n为样本量大小，k是预测变量数目
cutoff <- 4/(nrow(states)-length(fit$coefficients)-2)
plot(fit,which=4,cook.levels=cutoff)
abline(h=cutoff,lty=2,col="red")
image.png
  1. 变量添加图
    变量添加图是针对每个预测的变量Xk，绘制Xk在其它k-1个预测变量上回归惨纸质相对于相应变量在其它k-1个预测变量上的回归的残差值，在这个图中我们可以看到，如果改变或删除某个数值，如何影响线条的斜率，我们就知道这些点的影响力，如果一个点对直线极其强大，那么这个点属于强影响点。
library(car)
avPlots(fit,ask=FALSE,id.method="identity")
image.png

比如Income中的Alask这个点，删除这个点的话，直线应该是右旋转，斜率变成负值，这个点的影响力还是挺大的。

异常点（离群点、杠杆值和强影响点）的综合判断，通过influencePlot()函数

image.png

● 实际上还是通过标准差-2到2之间的数值来判断是否是离群点
● 帽子值超过0.2或0.3认为是高杠杆值点
● 强影响点应该是通过圆圈的大小，可以看到强影响的点应该是有两个Nevada和Alaska

2.4 回归模型的改善

根据上面的我们可以进行更好的回归诊断
实际上也是这些常见的，就是删除异常值的观测值，添加或者删除一些变量，变量转换以及使用其它的回归方法。

2.5 最佳的回归模型模型（模型比较和变量选择）

模型比较

  1. anova函数比较两个灿涛模型的拟合优度


    image.png

    可以看出来模型2嵌套在模型1中，但是最终的检测结果是F指是0.994，也就是不需要将这连个变量添加到模型中。

  2. AIC比较模型
    AIC比较了模型的统计拟合度和用来拟合的参数数目，AIC值比较小的模型要优先选择。


    image.png

变量的选择

实际上从很多候选变量中挑选最终的预测变量有两种方法：逐步回归法和全子集回归

  1. 逐步回归
    逐步回归，模型一次添加或删除一个变量，知道达到某个判定为准。
    逐步回归分为了三种：向前逐步回归、向后逐步回归和向前向后逐步回归。
    ● 向前逐步回归，每次添加一个预测变量到模型中，知道添加变量不会使得模型改善。
    ● 向后逐步回归，模型一开始包含所有的变量，逐一删除直到会降低模型质量
    ● 向前向后逐步回归，变量灭磁进入一个，每一步中，模型会重新评价，对模型没有贡献的变量删除，预测变量会反复的删除和提出，直到获得最佳模型。
library(MASS)
states<-as.data.frame(state.x77[,c("Murder","Population","Illiteracy",
                                   "Income","Frost")])
fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)
stepAIC(fit,direction="backward")
image.png

但你看，实际上，通过这种方式可能找到一个好模型，但是没有办法保证是最佳模型，因为你逐一删除，而不是所有模型都评估了，所以全子集回归法应运而生

  1. 全子集回归
    全子集回归检验所有可能的模型。
    通过leaps包中的regsubsets()函数实现，可以通过R平方、调整R平方或Mallows Cp统计量选择最佳模型
    ● R平方是预测变量解释想想变量的程度
    ● 调整R平方与其类似，但是考虑了模型的参数数目。R2可能丢失数据的偶然变异信息，调整R平方提供真实的R平方估计
    ● Mallows Cp统计量可作为逐步回归的判定标准，好的模型的Cp统计量接近模型参数数目(包括截距项)
    方法1：
library(leaps)
states<-as.data.frame(state.x77[,c("Murder","Population","Illiteracy",
                                   "Income","Frost")])
leaps <- regsubsets(Murder~Population+Illiteracy+Income+Frost,
                    data=states,nbest=4
                    )
plot(leaps,scale="adjr2")
image.png

此处，可以看到双预测变量模型（Population和Illiteracy）是最佳模型，其调整的R2最大，颜色最深，，代表了好的模型。

library(car)
subsets(leaps,statistic="cp",main="Cp Plot for All Subsets Regression")
abline(1,1,lty=2,col="red")
image.png

约好的模型，越是接近截距项和斜率为1的直线。
上面讲述了模型评估好坏和变量的选择，实际上哈哈，我个人认为这一张应该放在前面，就是说比如一个数据，我想通过模型进行解释，首先我要选择模型的类型，比如选择哪个模型进行尝试，泊松模型、Logistics模型，非线性模型等，还是其它的；上面的例子都是我们手动选择一些变量进行模型的建立，但实际上，我们应该先通过全子集回归模型或其他的模型然后来选择一定量的参数，选择好的模型，然后进行回归诊断，评估模型的统计假设是否成立，如果不行，剔除一些数据，做变量替换等等。OK，如果这个模型使用这个数据，为了更好的解释模型，或者其他情况，我们可以进行模型的泛化能力评估等。

2.6 深层次分析（模型的泛化能力和变量的相对重要性）

模型的泛化能力

实际上，我们通过模型是解释了该数据集，但是对于该模型能否解释其他的数据集或者进行预测，我们不得而知，这个时候我们就需要评估模型的泛化能力，其能否进行很好的预测。这时候用的是交叉验证。
所谓的交叉验证，就是取一定比例的数据作为训练样本，一部分样本作为保留样本，通过训练样本获得回归方程，通过保留样本进行预测。
k重交叉验证，就是分为k个子样本，轮流将k-1个子样本作为训练集，另一个子样本作为保留集，获得k个预测方程，记录k个保留样本的预测表现结果，与其平均机和。
bootstrap中的crossval函数可进行k重交叉验证

shrinkage <- function(fit, k=10){
require(bootstrap)
theta.fit <- function(x,y){lsfit(x,y)}
theta.predict <- function(fit,x){cbind(1,x)%*%fit$coef}
x <- fit$model[,2:ncol(fit$model)]
y <- fit$model[,1]
results <- crossval(x, y, theta.fit, theta.predict, ngroup=k)
r2 <- cor(y, fit$fitted.values)^2
r2cv <- cor(y, results$cv.fit)^2
cat("Original R-square =", r2, "\n")
cat(k, "Fold Cross-Validated R-square =", r2cv, "\n")
cat("Change =", r2-r2cv, "\n")
}
image.png

可以看到初始样品的R2为0.567，但是经过10重交叉验证，也就是泛化能力的R2是0.4481，但是我们最好挑选有很好泛化能力的模型。
比如重新验证


image.png

该模型的10重交叉验证的泛化能力0.53的R2，还是不错的，所以可以选择这个模型。

变量的相对重要性

有的时候研究者实际关心在一堆变量中那些变量是最重要的。我们可以通过一些方法来进行检测。

  1. 比较标准化的回归系数
    表示其它预测变量不变时，该预测变量一个标准差引起相应变量的预测变化
states<-as.data.frame(states.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
zstates<-as.data.frame(scale(states))
zfit <- lm(Murder~Population+Income+Illiteracy+Frost,data=zstates)
coef(zfit)
image.png

可以看到最重要的两个变量是Illiteracy和Population。

  1. 将相对重要性看做每个预测变量对R平方贡献
  2. relaimpo包包含了一些相对重要性的评价方法
    相对权重，是对所有可能子模型添加一个预测变量引起的R平方平均增加量的一个近似值。
relweights <- function(fit,...){
R <- cor(fit$model)
nvar <- ncol(R)
rxx <- R[2:nvar, 2:nvar]
rxy <- R[2:nvar, 1]
svd <- eigen(rxx)
evec <- svd$vectors
ev <- svd$values
delta <- diag(sqrt(ev))
lambda <- evec %*% delta %*% t(evec)
lambdasq <- lambda ^ 2
beta <- solve(lambda) %*% rxy
rsquare <- colSums(beta ^ 2)
rawwgt <- lambdasq %*% beta ^ 2
import <- (rawwgt / rsquare) * 100
import <- as.data.frame(import)
row.names(import) <- names(fit$model[2:nvar])
names(import) <- "Weights"
import <- import[order(import),1, drop=FALSE]
dotchart(import$Weights, labels=row.names(import),
xlab="% of R-Square", pch=19,
main="Relative Importance of Predictor Variables",
sub=paste("Total R-Square=", round(rsquare, digits=3)),
...)
return(import)
}

上面这个代码是我直接拷贝的，我没看。

> states <- as.data.frame(state.x77[,c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
> relweights(fit, col="blue")
Weights
Income 5.49
Population 14.72
Frost 20.79
Illiteracy 59.00
image.png

2.7 总结

这部分回归内容介绍了回归模型有哪些，回归模型的形式是什么，以线性模型举例说明（一元线性、多项式、多元线性等），回归诊断，异常值等、提高回归方程能力，模型的评估和变量选择，以及最终的模型的泛化和变量的相对重要性。
实际上，从我个人的角度出发，我认为，首先对于个人的一套数据集合而言，可能要通过以下的方式选择模型：

  1. 初步选择一个模型(根据数据的类型初步选择一个模型）
  2. 初步选择一个模型后，我们应该看一下变量中哪些最重要，我们可以通过比较标准化的回归系数或者是通过全子集回归或者逐步回归来挑选一些重要的变量。（而实际R语言实战这本书中这部分内容是放在后面的，模型吃的初建是通过认为选定变量挑选的）
  3. 然后根据挑选的变量来构建模型，构建模型后评估该模型的统计假设。当然了本文中模型是线性模型，要求符合正态性、独立性、线性、方差同质性等。如果不符合条件或者存在离群点，可以检测删除，重建模型。
  4. 模型构建完成后，模型重建后可再次进行回归诊断，评判统计假设。
  5. 而后判定模型的泛化能力，挑选泛化能力不错的模型（前提是模型能比较好的拟合数据）
    实际上，很多数据不满足统计假设，比如正态性，而在12章提供了置换随机化实验，可以不进行正态分布等。
©著作权归作者所有,转载或内容合作请联系作者
平台声明：文章内容（如有图片或视频亦包括在内）由作者上传并发布，文章内容仅代表作者本人观点，简书系信息发布平台，仅提供信息存储服务。
禁止转载，如需转载请通过简信或评论联系作者。

推荐阅读更多精彩内容