Video在线视频

下面我们开始为大家介绍一下怎么用R语言

来进行一个多元线性回归

当然又是一个很简单的介绍了

我们打开我准备好的

第十一周的新text这个新text你也可以从我们学堂在线

下载下来

我双击它在mac里面

直接就可以打开这个

R的平台了

那

用到的例子仍然是我上节课

就介绍过那本书的例子啊

这个数据叫做Boston

如果有兴趣的同学

欢迎大家下载那本书

然后到它的lab里面

那里面有专门有lab的介绍啊

来看看整个它具体的每一个细节的操作过程

那我们来继续看一下这个数据啊

这个数据仍然是想找到

怎么样可以解释

房屋价值的中位数

就是房子的均价

然后它这里面有13个predictor啊

比如说它这里面叫做rm average number of rooms per house

然后age average age of houses 房龄

然后lstat percent of households with low socioeconomic status

lstat这个因为上节课就有用到啊

它是

在这一栋房子里面

低社会经济地位住户的占比

那我们上节课就建立了一个

特别简单的线性回归

想看看这个

附近低社会经济地位人的比例和

和这个house value和房屋的价值有没有关系啊

那

重复一下上周的操作啊首先

因为数据是来自于叫MASS mass

这个程序包的

这个贡献包的

捐献包的首先我们要library

把这个捐献包把它叫出来

必须要把

它叫出来以后我们才可以

把这个数据读进来

然后后面我们还会遇到

ISLR这个捐献包的数据啊

所以我现在也先把它library把这个

捐献包叫出来

然后我attach

Boston这个数据

Boston在mass里面啊

我把它attach下来

后面你再

指定变量名的时候R就知道

到哪里去找了啊否则的话它就会报错

那我

上周就演示过的这个简单的线性回归用到的

是lm 这个function 啊

我说简称它就是linear model 线性模型的意思

然后前面这个lm.fit

这个是你自己可以随便任意取名的啊只要你自己

认识就可以我管它叫

a linear model拟合lm.fit它=lm

然后你的因变量y放在左边

是 median house value medv

然后它和谁形成关系呢

低社会经济地位人群的占比lstat

然后我现在

运行它一下它就默默的

这个跑下来了啊

因为线性回归在上节课我已经讲了我就

不再给大家看它的结果了

我们这节课想讲一下什么呢

就是做完任何一个线性回归以后其实

你可以用一个图

很快的进行一些诊断

我们刚刚讲了说在模型的

拟合过程中

可能会遇到很多的潜在问题

比如说我会担心呀我的xy之间是不是一个线性的关系呢

我还会担心

模型里面有没有异常值

或者是哪一个x它会对y产生一些不良的影响

那其实R语言呢

很方便的就可以给我们出来很多有用的图啊

那我下面这个命令

如果有经验的同学你们也可以认出来啊

认出来了我现在用parameter就是要画图了啊

par这个命令后面跟

mfrow

实际上我生成了一个

2*2的图形矩阵

就是我想

接下来我想画四个图

这个图

这四个图的排列方式呢就是两行两列

那我现在

回车一下你看它生成了一个

一页白纸啊实际上这是一个

2*2的白纸

然后下面我们用plot

lm.fit lm.fit就是我刚才

对回归这个对象操作的命名lm.fit

你看我只是plot让它画一个图那它就很自觉地画出了

4张非常fancy非常美丽的图啊

首先我们看左上角第一个图啊

这左上角第一个图是我们

之前讲这潜在问题的时候

解决问题的神器一样的那个

残差图

我们说如果

你的模型没什么问题的话

残差

应该随机的

散落在这个残差图的各处

然后如果我们对它进行一个比较实际情况的拟合

就中间那条红线啊做一个平滑的话

这条红线应该和0

距离很近

那你现在看到的这个图形呢

显然我们看残差没有

残差

并不是散落各处的啊感觉是有一个趋势的

而且

在你的拟合值是

20-25的时候它好像形成了一个聚集

所以不是很美丽

然后你看左下角这个图啊

它的y轴管它叫做standardized residuals

其实就是我们之前说的studendized residuals

就是

进行了某种标准化操作之后的残差图

仍然我也是希望这个残差图

它里面的点是散落各处的那

好在这个标准化以后的残差图看似

比直接的这个残差图稍微有了些进步

然后右上角是

standardized residuals的Q-Q图我们关注的是

我的残差项是不是服从正态分布

如果说残差项服从正态分布的话

我希望这些小圈圈应该都落在对角线上

就是我的theoretical quantiles 和我的standardized residuals

应该是一样大的

值是一样的啊而不是像我们现在看到这样一个尤其是在尾部

我们看这些点和这个虚线还是有一定距离的

那右下角residuals 和leverage

residuals就是残差

leverage是杠杆啊

我们这里面

显示的是cook's distance啊

然后仍然

纵轴是standardized residuals

关注的是有没有一些

点它会带来特别的杠杆效应

这是我特别简单的给大家讲解一下这个

模型的诊断啊

然后如果说我想一个图一个图的看也是可以的

我可以第一个你看plot

用plot这个命令呢

我可以y轴

对应的是 predict（lm.fit） 这是什么意思

我让它画的是

lm.fit这个回归里面的预测值predict

然后x是residuals

我现在运行一下

它实际上显示的就是左上角那个图

然后以此类推

有lm.fit然后rstudent

这是

studendized residuals

还有plot hatvalues

hatvalue大家还记得是什么吗

是看有没有

杠杆作用的那个h值

然后which.max

这个hatvalue那个最大它就有问题了啊

which.max你们

看不到图形了我们看到这个R的这个

这个console窗口它会告诉你说

第375个观测值的hatvalue是最大的

那这个时候你可以回到原始数据

看看第375个值

有没有什么可疑之处

这是

仍然是简单线性回归啊

我们看一看关于模型诊断的问题

然后下面我们就要进入到

多元线性回归的内容了啊

回归我们说

回归的时候并不总是只有一个

自变量经常情况下我们想关注

在控制了很多个自变量情况下某一个变量和

因变量的关系

那这里面呢我们

关于这个

这个例子的时候我们不说理论啊我们只是

讲解一下

R的一些功能所以我这是随便的比如说我把

age放进去age代表的是房龄

对吧我想看看房龄

对于整个

房屋的均价有没有影响

所以在原来的基础上我只是

加入了一个变量

加入一个变量操作非常简单直接在模型里面打

+age

就可以了然后

你也可以以防万一啊告诉

R说我要用到的数据是

Boston这个数据

lm.fit

然后我现在要summary给我展示

一下你的回归结果

用summary这个

summary这个函数展示回归结果你看

它会告诉你说

lstat它前面的系数是

-1.03

这说明

低社会经济地位的人占比越高

房屋价值就越低这个我们的

预期是差不多的而且它是一个非常显著的一个

影响因素或者它是

显著的和y值相关

然后age呢

说明好像

房龄越长

房屋的均价

是越高的

当然这个是两星显著啊

并没有显著性并没有

lstat这么好

那下面我刚才说我这个

数据里面一共有13个自变量啊

我特别懒

我直接要把这十三个

自变量一股脑的扔进模型里可不可以

仍然也是可以这样操作的啊

仍然用lm这个函数y不变啊

median value medv

然后这个破折号不变

后面加一个点这个点

经常代表是空值啊就

我什么都不敲

进去的情况下它会把

Boston这个数据里面除了

medv所有的变量都

被你一股脑的扔进模型里面

然后我summary一下你可以看看

是不是这样操作的

你们看这样犯罪率啊什么rm啊

disrad全都扔进去了

然后

我们还幸运地发现好像大多数的

自变量都和y

产生了某种显著的相关关系啊

然后下面呢我稍微再说一下关于

lm的这个对象啊

我们用lm这个函数进行多元回归

当然除了我们显示出来的

这些个输出值啊summary给了我们很多

模型回归的结果

但除了这些值以后实际上还有很多很多结果

它为了怕烦到你

它不全都输出出来了

那有的同学呢

可能你有一些特别需要需要报一些特别的

统计量总结性统计量

那你还需要在你的这个模型里再进一步挖掘一下

我这里面其实lm

给大家展示出来的东西一般是你做实际的研究

该报的这些结果都已经

囊括了

那如果说以后你做一些

更加高级的研究更加稀有的研究需要

报一些很特别的

总结性统计量的时候

你可以

首先啊用

用一个问号

然后summary.lm

就想知道

lm这个函数

除了summary以外

除summary产生这些output以外还有没有

其他的output啊

我现在打一个问号你可以看到

它会说

这个description啊它summary method for class“lm”

然后

这里面有一个value啊

它会告诉你说实际上在lm这个函数

这个函数被运行以后它产生了很多很多

统计量比如说包括residuals

残差

coefficients aliased

然后σ df

很多很多的这个总结性统计量

所以说有兴趣的话你们自己可以

好好的阅读一下

我举一个例子你看着里面

有σ

这个σ其实就是我们的rse啊

如果你看它的阅读

假设说

本来的这个summary默认的输出里面没有包括

σ这个数

你可以通过下面这个命令把这个结果叫出来

summary( lm.fit)

然后这边的啊

然后加一个美元$sigma

看

它给你

结果说

summary (lm.fit) sigma=4.74

它刚好和上面我们结果里面这个

residual standard error4.745是一模一样的啊

以此类推

比如说你想把r方

交出来呢就用summary （lm.fit） 然后$r.sq就可以了

然后下面还有一个

输出值有时候我们会关注啊叫做

膨胀因子

对吧如果你还记得

有膨胀因子的话

膨胀因子的话

代表的是模型里面有没有共线性

问题 它是这样一个指征啊

那并不是

就是R的标配里面

它并没有装膨胀因子的函数

所以为了找到这个

非常棒的函数呢我需要另外的在

安装一个捐献包叫做car

当然我已经安装好了然后

下面我要用library再把这个car叫出来

你们可能没有安装这个捐献包啊

这种情况下一定再用install.packages（）

把它先装进去

然后下面用到的函数

装进去了以后用到的函数就是

vif （lm.fit）我们看一下啊

直接就显示出了

每一个自变量它所对应的膨胀因子

然后

我们说如果说膨胀因子大于10了

那就说明

这个变量膨胀的有点厉害

共线性问题可能就有了那你就需要解决一下

然后我们这里面的变量呢

好像问题不是太大啊

有比如说

稍微有一点问题的 rad 还有tax 相对比较高

值得你关注一下

然后下面

再给大家介绍一个命令就是

如果说我刚才13个自变量一股脑地扔到这个

多元线性回归模型里面了

如果说里面有一两个变量

我不想要比较简单的操作

你可以直接在原来这个

曲线点的基础上

减去这个变量

我们看刚才这个模型

age这个变量在你加入了其他

变量以后它变得不显著了

你看我点亮的这一行啊

没有星星了

那如果我想把age

扔出去我直接用medv曲线点

减age就可以了

然后从新summary

就是一个没有age这个模型看起来

很厉害除了这个

indus这个变量以外其他全是显著的

然后这还有一种做法啊在原来的模型基础之上

敲lm.fit然后

破折号点减age也可以

好刚刚讲完了基本的多元回归的建模我们来

探讨一些小细节问题啊

假设说你的模型里面想加入交互项

那R也是可以很简单的操作的啊我们

前面很简略的

讲了一下什么是交互项啊

就是我有的时候想关注一下

两个自变量会不会

这个惺惺相惜的产生一些更

剧烈的在原有基础上

更强大的影响

那假设刚才这个例子啊

我的模型里假设

仍然只有lstst和

age这两个变量

现在我想

加入lstat和black

这个变量的交互

那交互很简单

只是在两个变量之间加入一个冒号

就可以了啊

我们假设这一个

变量我直接summary （lm())一气呵成

让它一下显示出来了啊

我们看这个结果

模型里面lstat

和black的交互项

它前面的系数呢等于0.0001474

然后P值是0.346发现这个交互项其实是不显著的

没有任何作用

然后下面仍然也是一种加入交互项的一种

方式这种方式

实际上一种更偷懒的方式

方法了啊

前面那个冒号只加了交互项

这个交互项本身那如果说

我偷懒我想一次性的把

自变量加进去

自变量以及它们的交互项一起加进去

我直接用星号的方法啊

lstat星号age

相乘的意思啊

我并没有打之前的lstat+age了

我在运行一下你看这时候lstat age

以及它们的交互项就一起扔到了模型里面

这个是怎么用R加交互项

然后下面

我们之前谈到如果说x和y

不是线性关系我就

需要对它进行一下线性转换

在

linear model这个函数的环境下

进行线性转换是一个特别简单的事情

那怎么做呢

假设说我还是和lstat这个

低社会经济地位人群占比

这个变量我就和它干上了

然后我要

增加一个变量我要

加入lstat的平方

那

有些同学可能直接想

在模型里面加lstst然后一个

小尖儿

然后乘以2

就好了就加平方的方式

但是在模型

在建模的这个环境之下呢

这个

小尖儿它还是有其他的作用的

所以你前面要加一个大写的

I把这个项目括起来

否则的话它会报错的啊

然后lm.fit2

等于linear model这个函数

加上后面那一坨东西然后我看看summary啊

summary告诉我们说lstat和

它的平方项它的二次项

都是特别显著的

那可能这就告诉我们说

做一下线性转换是有好处的

那下面我就想检验一下所谓的

非线性转换是不是有

是不是对原始的模型有了一个显著的作用

你可以

比较一下有它和没它的两种模型啊

没它的话我直接还是原来的

那个一元线性模型啊

下面有一个特别好用的命令可以比较

某两个模型之间有没有显著的差距

用到的命令是anova anova是方差分析了啊

这么anova如果把你的每个模型的这个

赋值放进去啊lm.fit

这是

一元线性模型

lm.fit2代表的是

有二次项的那个模型如果我想

比较这两个模型哪个好直接用anova

它的原假设你们猜猜是什么

原假设必然是我们说什么事都没发生

原来怎么样现在还是怎么样那

在anova这个检验里面

实际上原假设就是两个

模型的拟合程度一样好

加入了某一个新的变量

对于原来的模型没有显著的影响

所以原假设是没影响两个模型一样好

那你一定期盼的是

推翻原假设看到一个显著的P值

那我先在运行它一下

你会看到这里面F检验

它的p值是小于

2.2乘上10的-16次方

非常小那说明加入了

平方项以后

模型有非常显著的提高

然后下面

最后一个问题了啊最后一个小问题我们有一个虚拟变量的问题

就是假设说你的模型里面有

叫定性型变量

分类型变量你要怎么处理

那我们看

这时候我们用到

书里面啊

之前我们看library (ISLR)这个捐献包里有一个数据

叫做Carseats 就是

汽车

座位的销量这么一个例子

我们看 fix(Carseats) 这是我们数据啊

对于数据里面有很多的变量有

CompPrice Income Advertising Population Price

如果对这个数据

感到非常费解

我都不知道这数据说什么怎么办呢

讲了很多次了啊

我们做法是问号然后 Carseats

它会给你一个解释啊它说这是

A simulated data set

这是一个模拟的数据虚拟的数据

它包括了

sales of child car seats就是

儿童那个安全座椅啊

400个不同的商店销售儿童座椅

的情况

然后下面这个

dataframe这个数据里面呢有11个变量

然后它会给你解释每一个变量分别是什么

然后我可以用names这个函数

让它很快地告诉我说我这个数据里面有哪些变量它们分别叫什么名字

下面我要把这个数据attach进去啊

然后下面你在找变量的时候

不用再告诉R去哪一个数据找了啊

然后

假设说我需要找关系啊我需要

找什么关系呢我的

因变量是sales销售情况卖出去了多少

然后我的自变量

假设说

我想看

所有的自变量

以及还有一些交叉项的情况啊

我看直接用lm.fit =lm

sales因变量然后一个小曲线 点

点的话就相当于

把刚才的数据中

除了sales所有的变量都一股脑的扔进去了

同时我还加入了两个交互项

一个是income和广告的交互项一个是price和age的交互项对吧

我lm.fit直接就跑出来了然后summary一下

然后summary一下得到了一个结果

然后你关注一下我这里面有一个

变量叫做

shelvelocation

shel

veloc

然后这里面是good然后medium

这个变量我们再回到看看原来数据啊

我们看原来数据这里面

这一列我点亮了这一列

就是

它代表的是这个

安全座椅在这个店里面

摆放的位置好不好或者是占的地方大不大

如果安全座椅在

摆在店里非常鲜明的位置占的

占地很大的话

人们很快就注意到它想看看它会不会影响销量啊

然后这是一个

分类型变量我并没有一个很具体的数值

我只是把

这个空间这个

货架好坏分成了

bad medium 还有good

差中好三级

这是一个典型的分类变量啊

那放到模型里显然要把它当虚拟变量处理

对于一个有三种分类取值的这样一个

分类变量我们

会生成两个新的变量

你想想两个新的变量

假设说对于一个bad good medium

我其实什么都没有操作而自动已经把它

当成分类变量处理了

而且它非常自觉地给你建了两个变量

分别是good和medium

然后参照组显然就是那个

bad了啊就是坏是参照组

那这个good什么意思呢

这是坏的那组和好的那一组

相比啊如果你把这个

安全座椅放到非常好的销售位置的话

哇它销量有特别大的增加啊

有多不同呢它们的不同是4.84

个单位

然后 medium和差之间的销量不同在

1.95个单位啊这是怎么去解读它啊我们说

这个4.84代表的是

good和bad之间销量的差距

销量的差距啊那

我现在想看看你到底这个虚拟变量怎么

构建的我可以用一个函数叫做contrasts

contrasts()

contrast good medium

原来的变量的取值是bad good

然后在

建立虚拟过程中我就剩下

两个变量了啊

这怎么理解呢

假设说原来的变量

分类是在bad是在不好的货架时候

新的虚拟变量这个

good这个变量的取值就是0

然后 medium的取值也是0

然后如果原来是good那显然

good这个取值就是1medium是0

然后 medium的话

如果本身是medium

那么新变量里面good取值就是0medium就是1

所以这就告诉你们其实

两个新建的虚拟变量

足够可以代表一个原来分类是

三分类的这么一个变量了

那

实际上为什么R可以

这么自觉的就把它处置成分类变量呢

因为我这个数据输入的时候它的模式就是factor模式

我可以用class这个命令来看看我所有变量的

所有变量的类型

我刚才有那么多个变量啊假设说我这里面有age那个变量

它会告诉你说age这个变量是数值型的

以此类推那我试试

shelveloc

它会告诉你说这是factor分类型变量

所以如果说你的变量

是分类型变量你可以用这个命令假设说

你输入时候数字

并不是我现在这样输bad good medium 用真的

字符来表示的

假设说你的bad是1good是2 medium是3

那R可能没有那么聪明

他就会以为你这是个数值型变量这个时候怎么办

你要用as.factor()

这个命令把它改过来啊

比如说如果说

ShelveLoc刚好是

数值型的输入形式

那你需要用 as.factor(ShelveLoc)把它改过来

当然这是数据清理啊数据准备的时候

非常细节的操作好多好多的事情你要做

如果说你手头正面临着特别大的一套数据啊

一定要

清晰地知道你输入的时候每一个数据的

存储类型它到底是数值型

分类型还是什么其他的类型啊

好

这是我们没有多少时间啊这是非常非常

简略的给大家介绍一下R里面在

进行多元回归的一些时候的

特别基本的功能

如果你有兴趣

做一些更加细节的操作

一个是看我刚才推荐那本书

另外一个多打一些问号

对于我这里面用到的每一个

函数你都可以打一个问号

看一看具体还有哪些更加

精细的功能

好这就是本周的R语言的演示