多元线性回归
- 确定因变量y和影响因变量的k个自变量
假定为线性关系
,建立线性关系模型- 对模型进行估计与检验
- 判别是否有
多重共线性
,如果存在,进行处理 - 利用回归方程进行预测
- 对回归模型进行诊断
回归模型拟合
mode<-lm(y~x1+x2+x3+x4+x5,data=table)
summary(mode)
计算回归系数的置信区间
confint(mode,level=0.95)
输出方差分析表
anova(mode)
绘制残差图诊断模型
plot(mode)
对模型中各对自变量之间的相关系数进行检验
library(psych)
corr.test(table[3:7],use = "complete")
容忍度与膨胀因子
library(car)
vif(mode) #膨胀因子
1/vif(mode) #容忍度
逐步回归
mode1<-step(mode)
标准化回归系数
library(lm.beta)
mode.beta<-lm.beta(mode)
summary(mode.beta)
anova(要求模型相互嵌套)
model1<-lm(y~x1+x2+x3+x4+x5,data=table)
model2<-lm(y~x1+x2+x5,data=table)
anova(model1,model2)
AIC(不要求嵌套)
model1<-lm(y~x1+x2+x3+x4+x5,data=table)
model2<-lm(y~x1+x2+x5,data=table)
AIC(model1,model2)
文章目录
10.1 多元线性回归模型及参数估计
10.1.1 回归模型与回归方程
- 多元线性回归模型
e为误差项
对e的三个基本假定
- 正态性
- 方差齐性
- 独立性
- 多元线性回归方程
- 估计的多元线性回归方程
10.1.2 参数的最小二乘估计
- 回归模型拟合
> table<-read.csv("/Users/zhourui/Documents/example10_1.csv")
> mode<-lm(y~x1+x2+x3+x4+x5,data=table)
> summary(mode)
Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x5, data = table)
Residuals:
Min 1Q Median 3Q Max
-16.7204 -6.0600 0.7152 3.2144 21.4805
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.2604768 10.4679833 0.407 0.68856
x1 0.1273254 0.0959790 1.327 0.20037
x2 0.1605660 0.0556834 2.884 0.00952 **
x3 0.0007636 0.0013556 0.563 0.57982
x4 -0.3331990 0.3986248 -0.836 0.41362
x5 -0.5746462 0.3087506 -1.861 0.07826 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 10.65 on 19 degrees of freedom
Multiple R-squared: 0.8518, Adjusted R-squared: 0.8128
F-statistic: 21.84 on 5 and 19 DF, p-value: 2.835e-07
y = 0.1273254 * x1 + 0.1605660 * x2+0.0007636 * x3 -0.3331990 * x4 -0.5746462 * x5+4.2604768
x2对y值影响较大
- 计算回归系数的置信区间
> confint(mode,level=0.95)
2.5 % 97.5 %
(Intercept) -17.649264072 26.170217667
x1 -0.073561002 0.328211809
x2 0.044019355 0.277112598
x3 -0.002073719 0.003600932
x4 -1.167530271 0.501132297
x5 -1.220868586 0.071576251
- 输出方差分析表
> anova(mode)
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x1 1 10508.9 10508.9 92.7389 9.625e-09 ***
x2 1 1347.1 1347.1 11.8878 0.002696 **
x3 1 85.4 85.4 0.7539 0.396074
x4 1 40.5 40.5 0.3573 0.557082
x5 1 392.5 392.5 3.4641 0.078262 .
Residuals 19 2153.0 113.3
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
x1、x2对y值影响较大
10.2 拟合优度和显著性检验
10.2.1 模型的拟合优度
用多重决定系数
评估
1. 多重决定系数
R方
Multiple R-squared: 0.8518
调整后的R方
考虑样本量
和模型中自变量的个数
Adjusted R-squared: 0.8128
2. 估计标准误
Residual standard error: 10.65 on 19 degrees of freedom
10.2.2 模型的显著性检验
- 在一元线性回归中
只有一个自变量,F检验(线性关系检验)
和t检验(回归系数检验)
是等价的
- 多元线性回归
F检验:检验因变量同多个自变量的整体线性关系是否显著(有一个就显著)
t检验:对每个回归系数分别检验,以判断每个自变量对因变量的影响是否显著
1.线性关系检验
2.回归系数检验
3.模型诊断
- 绘制残差图诊断模型
> par(mfrow=c(1,2),mai=c(0.8,0.8,0.4,0.1),cex=0.8,cex.main=0.7)
> plot(mode,which=1:2)
10.3 多重共线性及其处理
10.3.1 多重共线性及其识别
- 多重共线性:当回归模型中两个或两个以上的自变量
彼此相关
1. 多重共线性所产生的问题
- 变量高度相关时,可能会造成回归结果的混乱
- 对参数估计值的正负号产生影响
2. 多重共线性的识别和处理
- 对模型中各对自变量之间的相关系数进行检验
> library(psych)
> corr.test(table[3:7],use = "complete")
Call:corr.test(x = table[3:7], use = "complete")
Correlation matrix
x1 x2 x3 x4 x5
x1 1.00 0.74 0.88 -0.62 -0.28
x2 0.74 1.00 0.55 -0.54 -0.32
x3 0.88 0.55 1.00 -0.52 -0.29
x4 -0.62 -0.54 -0.52 1.00 0.10
x5 -0.28 -0.32 -0.29 0.10 1.00
Sample Size
[1] 25
Probability values (Entries above the diagonal are adjusted for multiple tests.)
x1 x2 x3 x4 x5
x1 0.00 0.00 0.00 0.01 0.47
x2 0.00 0.00 0.03 0.03 0.46
x3 0.00 0.00 0.00 0.04 0.47
x4 0.00 0.01 0.01 0.00 0.65
x5 0.18 0.12 0.16 0.65 0.00
To see confidence intervals of the correlations, print with the short=FALSE option
x5对x1 0.47
x5对x2 0.46
x5对x3 0.47
x5对x4 0.65
均大于0.05(P值)
x5与其他自变量间关系不显著
use = “complete”(一对多,形成矩阵)
use = “pair”(一对一)
- 考察各回归系数的显著性
当模型的
F检验显著
多个回归系数的t检验却不显著
(同10_1)
- 分析回归系数的正负号
回归系数正负号与预期相反
- 用容忍度和方差膨胀因子来识别
容忍度小于0.1,存在严重共线性
膨胀因子(VIF)越大,共线性越严重,VIF大于10,存在严重共线性
> library(car)
> vif(mode) #膨胀因子
x1 x2 x3 x4 x5
8.233159 2.629940 5.184365 1.702361 1.174053
> 1/vif(mode) #容忍度
x1 x2 x3 x4 x5
0.1214601 0.3802368 0.1928877 0.5874195 0.8517500
回归模型的共线性并不严重
10.3.2 变量选择与逐步回归
在建立模型前有选择的确定进入模型的自变量
,可以避免多重共线性的问题
1. 向前选择
- 引入F最大(P最小)
- 加入的变量一定会被保留
2.向后剔除
- 剔除F最小(P最大)
- 剔除的变量不会再进入
3. 逐步回归
-
- 上述两种方法相结合
-
- 变量有可能被踢出,也有可能重新加入
- 变量有可能被踢出,也有可能重新加入
-
AIC
- 逐步回归变量选择
> mode1<-step(mode)
Start: AIC=123.39
y ~ x1 + x2 + x3 + x4 + x5
Df Sum of Sq RSS AIC
- x3 1 35.96 2189.0 121.81
- x4 1 79.17 2232.2 122.30
<none> 2153.0 123.39
- x1 1 199.42 2352.4 123.61
- x5 1 392.54 2545.6 125.58
- x2 1 942.22 3095.2 130.47
Step: AIC=121.81
y ~ x1 + x2 + x4 + x5
Df Sum of Sq RSS AIC
- x4 1 78.22 2267.2 120.69
<none> 2189.0 121.81
- x5 1 445.69 2634.7 124.44
- x2 1 925.88 3114.9 128.63
- x1 1 1133.27 3322.3 130.24
Step: AIC=120.69
y ~ x1 + x2 + x5
Df Sum of Sq RSS AIC
<none> 2267.2 120.69
- x5 1 404.28 2671.5 122.79
- x2 1 1050.90 3318.1 128.21
- x1 1 1661.83 3929.0 132.43
在最后一步中去掉三个变量中的任何一个都会增大AIC,所以运行停止
重新拟合后进行分析:
- 拟合逐步回归模型
> mode12<-lm(y~x1+x2+x5,data=table)
> summary(mode12)
Call:
lm(formula = y ~ x1 + x2 + x5, data = table)
Residuals:
Min 1Q Median 3Q Max
-14.027 -5.361 -1.560 2.304 23.001
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.68928 6.25242 -0.270 0.78966
x1 0.19022 0.04848 3.923 0.00078 ***
x2 0.15763 0.05052 3.120 0.00518 **
x5 -0.56979 0.29445 -1.935 0.06656 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 10.39 on 21 degrees of freedom
Multiple R-squared: 0.8439, Adjusted R-squared: 0.8216
F-statistic: 37.85 on 3 and 21 DF, p-value: 1.187e-08
- 逐步回归的方差分析表
> anova(mode12)
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x1 1 10508.9 10508.9 97.3392 2.452e-09 ***
x2 1 1347.1 1347.1 12.4775 0.001976 **
x5 1 404.3 404.3 3.7447 0.066558 .
Residuals 21 2267.2 108.0
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
y = -1.68928 + 0.19022 * x1 + 0.15763 * x2 -0.56979 * x5
- 模型诊断
> par(mfrow=c(1,2),mai=c(0.8,0.8,0.4,0.1),cex=0.8,cex.main=0.7)
> plot(mode12,which=1:2)
10.4 相对重要性和模型比较
哪些自变量对预测更重要
10.4.1 自变量的相对重要性
评估自变量相对重要性:比较标准回归系数
标准化回归系数越大,该自变量对因变量预测越重要
- 标准化回归系数
> library(lm.beta)
> mode.beta<-lm.beta(mode)
> summary(mode.beta)
Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x5, data = table)
Residuals:
Min 1Q Median 3Q Max
-16.7204 -6.0600 0.7152 3.2144 21.4805
Coefficients:
Estimate Standardized Std. Error t value Pr(>|t|)
(Intercept) 4.2604768 0.0000000 10.4679833 0.407 0.68856
x1 0.1273254 0.3361822 0.0959790 1.327 0.20037
x2 0.1605660 0.4130034 0.0556834 2.884 0.00952 **
x3 0.0007636 0.1132753 0.0013556 0.563 0.57982
x4 -0.3331990 -0.0963203 0.3986248 -0.836 0.41362
x5 -0.5746462 -0.1781104 0.3087506 -1.861 0.07826 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 10.65 on 19 degrees of freedom
Multiple R-squared: 0.8518, Adjusted R-squared: 0.8128
F-statistic: 21.84 on 5 and 19 DF, p-value: 2.835e-07
10.4.2 模型比较
简约模型:最佳模型
逐步回归方法虽然剔除不必要自变量,但是无法确认是否更好拟合,需要进一步比较
目的:找简约模型
- 嵌套模型
- 检验
- anova(要求模型相互嵌套)
> model1<-lm(y~x1+x2+x3+x4+x5,data=table)
> model2<-lm(y~x1+x2+x5,data=table)
> anova(model1,model2)
Analysis of Variance Table
Model 1: y ~ x1 + x2 + x3 + x4 + x5
Model 2: y ~ x1 + x2 + x5
Res.Df RSS Df Sum of Sq F Pr(>F)
1 19 2153.0
2 21 2267.2 -2 -114.17 0.5038 0.6121
RSS 残差平方和
Sum of Sq 交互效应残差和
P=0.6121,不拒绝H0,没有证据显示两个模型有显著差异
由于简约原则选择model2
- AIC(不要求嵌套)
> model1<-lm(y~x1+x2+x3+x4+x5,data=table)
> model2<-lm(y~x1+x2+x5,data=table)
> AIC(model1,model2)
df AIC
model1 7 196.3408
model2 5 193.6325
AIC 小说明模型用较少参数
逐步回归模型更好
10.5 利用回归方程进行预测
> pre<-predict(model2)
> res<-residuals(model2)
> zre<-rstandard(model2)
> x<-table[,c(3,4,7)]
> con_int<-predict(model2,x,interval = "confidence",level=0.95)
> pre_int<-predict(model2,x,interval = "prediction",level=0.95)
> musummary<-data.frame(营业额=table$y,点预测值=pre,残差=res,标准化残差=zre,置信下限=con_int[,2],置信上限=con_int[,3],预测下限=pre_int[,2],预测上限=pre_int[,3])
> round(musummary,3)
> round(musummary,3)
营业额 点预测值 残差 标准化残差 置信下限 置信上限 预测下限 预测上限
1 53.2 52.189 1.011 0.102 45.457 58.921 29.557 74.822
2 18.5 -4.501 23.001 2.359 -11.969 2.967 -27.363 18.361
3 11.3 21.963 -10.663 -1.072 15.685 28.240 -0.539 44.464
4 84.7 65.114 19.586 2.766 49.300 80.928 38.337 91.891
5 7.3 6.128 1.172 0.122 -2.283 14.540 -17.059 29.316
6 17.9 22.408 -4.508 -0.466 14.581 30.235 -0.574 45.390
7 2.5 1.278 1.222 0.125 -6.113 8.670 -21.559 24.116
8 27.3 34.719 -7.419 -0.747 28.400 41.038 12.206 57.232
9 5.9 10.628 -4.728 -0.472 4.904 16.351 -11.726 32.981
10 23.9 37.843 -13.943 -1.390 32.193 43.493 15.509 60.178
11 69.4 62.852 6.548 0.709 52.939 72.766 39.079 86.626
12 20.6 18.296 2.304 0.229 13.036 23.556 -3.943 40.535
13 1.9 -5.510 7.410 0.771 -13.694 2.674 -28.616 17.596
14 3.0 14.956 -11.956 -1.202 8.687 21.224 -7.543 37.455
15 7.3 8.860 -1.560 -0.169 -1.050 18.770 -14.913 32.632
16 46.2 29.831 16.369 1.699 21.745 37.917 6.759 52.903
17 78.8 78.016 0.784 0.112 62.105 93.927 51.182 104.850
18 11.1 13.167 -2.067 -0.209 6.420 19.915 -9.470 35.805
19 8.6 15.847 -7.247 -0.815 4.670 27.024 -8.481 40.175
20 48.9 50.727 -1.827 -0.184 44.205 57.250 28.156 73.299
21 22.1 27.461 -5.361 -0.536 21.671 33.251 5.090 49.831
22 11.1 0.233 10.867 1.354 -13.486 13.953 -25.362 25.829
23 8.6 22.627 -14.027 -1.395 17.181 28.073 0.344 44.911
24 48.9 48.961 -0.061 -0.006 42.722 55.200 26.470 71.452
25 22.1 27.005 -4.905 -0.490 21.220 32.790 4.636 49.374
- 求x为具体值时,日均营业额的点预测值、置信区间和预测区间(新值预测)
> x0<-data.frame(x1=50,x2=100,x5=10)
> predict(model2,newdata = x0)
1
17.88685
> predict(model2,x0,interval = "confidence",level=0.95)
fit lwr upr
1 17.88685 10.98784 24.78585
> predict(model2,x0,interval = "prediction",level=0.95)
fit lwr upr
1 17.88685 -4.795935 40.56963
10.6 哑变量回归
哑变量 == 分类类型变量
10.6.1 在模型中引入哑变量
10.6.2 含有一个哑变量的回归
> table<-read.csv("/Users/zhourui/Documents/example10_72.csv")
> model<-lm(y~X1,data=table)
> summary(model)
Call:
lm(formula = y ~ X1, data = table)
Residuals:
Min 1Q Median 3Q Max
-19.7604 -10.7832 0.7195 4.3343 28.9301
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.75023 5.25068 -1.095 0.285
X1 0.32394 0.04482 7.227 2.34e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 13.9 on 23 degrees of freedom
Multiple R-squared: 0.6943, Adjusted R-squared: 0.681
F-statistic: 52.23 on 1 and 23 DF, p-value: 2.343e-07
- 二元回归
> model<-lm(y~X1+x2,data=table)
> summary(model)
Call:
lm(formula = y ~ X1 + x2, data = table)
Residuals:
Min 1Q Median 3Q Max
-19.443 -11.579 -1.256 8.607 23.456
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.45413 4.69817 -1.799 0.08568 .
X1 0.28641 0.04145 6.909 6.15e-07 ***
x2方便 14.62088 5.17802 2.824 0.00989 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 12.17 on 22 degrees of freedom
Multiple R-squared: 0.7756, Adjusted R-squared: 0.7552
F-statistic: 38.02 on 2 and 22 DF, p-value: 7.269e-08
- 方差分析表
> anova(model)
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
X1 1 10085.9 10085.9 68.062 3.529e-08 ***
x2 1 1181.5 1181.5 7.973 0.009889 **
Residuals 22 3260.1 148.2
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
- 分析