学习笔记 | 回归模型 | 05 多变量回归分析
Multi-variable regression analyses
如果现在我需要做一个薄荷糖和肺功能直接关系的研究，你可能会产生以下怀疑
01 很可能你会认为：吸烟者更需要薄荷糖，吸烟者肺功能较弱，这很可能是主因
02 什么说服你这样去想，你会说，如果说不吸烟的薄荷糖使用者的肺功能比不吸烟不吃薄荷糖的人低，同理，吸烟的薄荷糖使用者的肺功能也比吸烟不吃薄荷糖的人低，我会更容易相信你
换句话说，为了证明，我需要提出吸烟状态再去考虑
有时候数据量的规模会导致无法使用简单的线性回归处理，那么，要怎么去对多个回归量进行预测？多个回归量又会导致什么？
01 有些回归量与Y值不相关
02 有些缺省的参数
这时我们就需要应用多变量的回归分析方法，其实就是先把数据进行分类再进行回归模型构建的方法。
练习1：
01 使用杀虫水数据（InsectSprays）数据概况：
> dim(InsectSprays)
[1] 72 2
> head(InsectSprays,15)
count spray
1 10 A
2 7 A
3 20 A
4 14 A
5 14 A
6 12 A
7 10 A
8 23 A
9 17 A
10 20 A
11 14 A
12 13 A
13 11 B
14 17 B
15 21 B
> summary(InsectSprays[,2])
A B C D E F
12 12 12 12 12 12
> sapply(InsectSprays,class)
count spray
"numeric" "factor"
02 生成count与spray的线性模型，R默认使用第一组作为截距参考，然后其它组则以此为参照计算t-test
这样会产生一个问题，就是Estimate列的数值并不是每一组的均值
> fit <- lm(count ~ spray, InsectSprays)
> summary(fit)$coef
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.5000000 1.132156 12.8074279 1.470512e-19
sprayB 0.8333333 1.601110 0.5204724 6.044761e-01
sprayC -12.4166667 1.601110 -7.7550382 7.266893e-11
sprayD -9.5833333 1.601110 -5.9854322 9.816910e-08
sprayE -11.0000000 1.601110 -6.8702352 2.753922e-09
sprayF 2.1666667 1.601110 1.3532281 1.805998e-01
于是用“-1”的方法处理
> nfit <- lm(count ~ spray -1, InsectSprays)
> summary(nfit)$coef
Estimate Std. Error t value Pr(>|t|)
sprayA 14.500000 1.132156 12.807428 1.470512e-19
sprayB 15.333333 1.132156 13.543487 1.001994e-20
sprayC 2.083333 1.132156 1.840148 7.024334e-02
sprayD 4.916667 1.132156 4.342749 4.953047e-05
sprayE 3.500000 1.132156 3.091448 2.916794e-03
sprayF 16.666667 1.132156 14.721181 1.573471e-22
练习2：
这个练习使用了WHO的饥荒数据，这些数据展示了全球范围的饥荒程度。
先来看一下这个数据hunger的一些基本情况：
> dim(hunger)
[1] 948 13
948个样例（行），13个指标（列）
> names(hunger)
[1] "X" "Indicator" "Data.Source" "PUBLISH.STATES" "Year" "WHO.region"
[7] "Country" "Sex" "Display.Value" "Numeric" "Low" "High"
[13] "Comments"
列Numeric的含义是有多少5岁以下的孩子的体重是低于标准的
下面对Numeric和Year生成回归模型
> fit<-lm(Numeric ~ Year, hunger)
> summary(fit)$coef
Estimate Std. Error t value Pr(>|t|)
(Intercept) 634.479660 121.1445995 5.237375 2.007699e-07
Year -0.308397 0.0605292 -5.095012 4.209412e-07
相关系数是 -0.30840，负相关表明随着时间，饥荒率在下降。
模型的斜率（634.479660）表示饥荒儿童在Year=0时的数值
下面我们来把男性和女性的数据过滤出来生成回归模型
> lmF <- lm(Numeric[Sex=="Female"] ~ Year[Sex=="Female"],hunger)
> lmM <- lm(Numeric[Sex=="Male"] ~ Year[Sex=="Male"],hunger)
通过图形化我们发现两条直线并不平行。
因此，我们用（性别+Year）来作为横坐标，如下
> lmBoth <- lm(Numeric ~ Year + Sex, hunger)
> summary(lmBoth)
Call:
lm(formula = Numeric ~ Year + Sex, data = hunger)
Residuals:
Min 1Q Median 3Q Max
-25.472 -11.297 -1.848 7.058 45.990
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 633.5283 120.8950 5.240 1.98e-07 ***
Year -0.3084 0.0604 -5.106 3.99e-07 ***
SexMale 1.9027 0.8576 2.219 0.0267 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 13.2 on 945 degrees of freedom
Multiple R-squared: 0.03175, Adjusted R-squared: 0.0297
F-statistic: 15.49 on 2 and 945 DF, p-value: 2.392e-07
SexFemale被当作是Intercept了，下面是Year和SexMale相对SexFemale给出的数值。
现在再把lmBoth图形化的话可以看到两条直线是平行的，因为他们有一样的斜率。
> lmInter <- lm(Numeric ~ Year + Sex + Sex*Year, hunger)| You are amazing!
> summary(lmInter)
Call:
lm(formula = Numeric ~ Year + Sex + Sex * Year, data = hunger)
Residuals:
Min 1Q Median 3Q Max
-25.913 -11.248 -1.853 7.087 46.146
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 603.50580 171.05519 3.528 0.000439 ***
Year -0.29340 0.08547 -3.433 0.000623 ***
SexMale 61.94772 241.90858 0.256 0.797946
Year:SexMale -0.03000 0.12087 -0.248 0.804022
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 13.21 on 944 degrees of freedom
Multiple R-squared: 0.03181, Adjusted R-squared: 0.02874
F-statistic: 10.34 on 3 and 944 DF, p-value: 1.064e-06
最后我们发现，当两个预测变量是连续的时候，处理起来会有点麻烦。
假设两个预测变量的其中一个是常量，那么模型的变化值就恒定了，解析如下：
当线性模型：Hi = b0 + (b1Ii) + (b2Yi)+ (b3IiYi) + ei
则 当I=5时，Hi就相当于一个等比数列，这时H单位变化量=b2+b3*5
