毕业设计实用模型（四）——回归模型的实现（R语言）

在毕业实用模型一¹、二²、三³中介绍了时间序列模型的建模思路与在R语言中的调参&＃xff0c;今天来讲解一下回归模型的建模思路和实现。
在这里你将会学到&＃xff1a;

• 推荐的回归模型的书籍
• 高效简介的相关性分析图
• 多元回归模型的建立
• 回归模型的诊断的阅读导读
• 有效变量的选择方法
• 回归常用的函数

本文的案例偏向R语言的实现和分析理论较少。在做回归模型之前&＃xff0c;介绍几本做回归的书籍&＃xff0c;其中一本是偏向理论的&＃xff0c;大学本科的教材&＃xff0c;有概率论数理统计基础即可看懂。一本是R语言实战的第八章&＃xff0c;偏重讲回归模型的实现。

1.1理论的书

要想看懂本文的思路&＃xff0c;需要一些基本的回归分析的理论&＃xff0c;下面这本书⁴可以给你这些&＃xff1a;考虑版权这里就不贴链接了。
还有一本书叫统计学习基础⁵&＃xff0c;这本书是偏向理论的统计模型&＃xff0c;介绍回归和相关的模型理论比较深一点&＃xff0c;如果遇到的问题上面的书解决不了可以选择求助这本书,本书的三、四、七、八章是关于回归的&＃xff0c;大家查阅学习。

1.2 R语言的书

R语言实战⁶被很多人说成R语言的圣经&＃xff0c;好好读完&＃xff0c;你的R语言水平会提高一个档次。今天主要参考的是他的第八章。

展示部分目录&＃xff1a;

做回归分析需要有实际的数据&＃xff0c;本篇使用随机模型的方式去生成规范化的数据进行建模。提前设定好真实参数&＃xff0c;可以更好的看到模型的效果。下面是构造数据的代码&＃xff1a;

library(MASS) # 生成多元正态分布
> n &＃61; 300; p &＃61; 6
> Beta <- c(rep(1,3),0,0,1) # 设置真实回归模型参数
> set.seed(0) # 设置种子
> x &＃61; mvrnorm(n &＃61; n, mu &＃61; rep(0, p-1), Sigma &＃61; diag(rep(1, p-1))) # 需要载入MASS包
> x6 <- x[,1] &＃43; rnorm(n, 0, 0.3)
> x &＃61; cbind(x, x6)
> colnames(x) <- paste0("x", 1:6)
> y &＃61; x %*% Beta &＃43; &＃43; rnorm(n, 0, 2)
> Data <- data.frame(x, y &＃61; y)
> head(Data)x1 x2 x3 x4 x5 x6
1 -0.9593164 -0.59188422 -1.6878010 -0.1244350 1.2629543 -0.9356879
2 -1.6203166 -0.37099306 0.6476460 1.4667446 -0.3262334 -1.6357428
3 0.8225133 0.08792426 0.4487942 0.6739287 1.3297993 0.7442811
4 0.1087127 -0.03472634 1.0263022 1.9564253 1.2724293 0.5785931
5 0.7609948 1.80637427 1.0749782 -0.2690410 0.4146414 0.6493902
6 -2.3062566 -0.34023607 0.4583096 -1.2445515 -1.5399500 -1.7836732y
1 -4.5420488
2 -2.5045986
3 0.9620528
4 3.0366120
5 1.3610465
6 -6.2472318


介绍一下数据的生成。六个自变量。其中$x6$是由$x1$线性生成的。也就是$x1、x6$具有很大的相关性。参数设置时。$x1、x2、x3、x6$的参数是1&＃xff0c;$x4、x5$的参数是0。一元回归的建立与可视化参见R语言可视化——ggplot2画回归曲线一文&＃xff0c;里面包含了一元回归的方差分析、回归线、R2R^2R2等R语言的实现。下面讨论多元回归模型的建立。

在相关性分析时&＃xff0c;往往需要绘制一幅比较清晰明了的相关图。下面我就对上述数据画出一部分吧相关图。

library(corrplot) # 画相关图的包
library(customLayout) # 加载拼图包
# 创建拼图画布
lay1 <- lay_new( mat &＃61; matrix(1:2, ncol &＃61; 1), heights &＃61; c(1,1)
)
lay2 <- lay_new( mat &＃61; matrix(1:4, ncol &＃61; 1), heights &＃61; c(1,1,1,1)
)
lay3 <- lay_bind_col(lay1, lay2, widths &＃61; c(2, 1)) # 合并画布
par(mar &＃61; c(1, 1, 1, 1))
lay_show(lay3)
lay_set(lay3)
M <- cor(Data)
col1 <- colorRampPalette(c("#7F0000", "red", "#FF7F00", "yellow", "white","cyan", "#007FFF", "blue","#00007F"))
corrplot(M, method &＃61; "number", col &＃61; "black", cl.pos &＃61; "n")
corrplot(M, method &＃61; "pie", order &＃61; "AOE")
corrplot(M, method &＃61; "number")
corrplot(M, order &＃61; "AOE", addCoef.col &＃61; "grey")
corrplot(M, order &＃61; "AOE", col &＃61; col1(20), cl.length &＃61; 21, addCoef.col &＃61; "grey")
corrplot(M, method &＃61; "color", col &＃61; col1(20), cl.length &＃61; 21, order &＃61; "AOE",addCoef.col &＃61; "grey")
plot(Data)


下面是两幅图的效果&＃xff1a;

可以在上面的相关图中按自己的需求选一幅进行分析。除此之外在建模中之前还应该&＃xff1a;

• 对$y$进行正态性检验。
• 对数据进行中心或者标准化&＃xff0c;尽量各个变量之间的的差距不是很大。具体理论参考文献⁴第75-78页。
• n是不是远大于p&＃xff0c;如果n不是远大于p或者数量比较接近&＃xff0c;则可能导致过拟合、估计的参数可信度不大。

4.1建立模型

> par(mfrow &＃61; c(2, 2))
> fit <- lm(y~., data &＃61; Data)
> summary(fit)
Call:
lm(formula &＃61; y ~ ., data &＃61; Data)
Residuals:Min 1Q Median 3Q Max
-5.8223 -1.3649 -0.1257 1.5422 6.0490
Coefficients:Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.04555 0.12520 -0.364 0.71627
x1 1.13727 0.40586 2.802 0.00541 **
x2 1.09968 0.11920 9.225 <2e-16 ***
x3 0.99557 0.12964 7.679 2.41e-13 ***
x4 -0.08748 0.12495 -0.700 0.48442
x5 0.03488 0.13029 0.268 0.78913
x6 1.02470 0.39657 2.584 0.01025 *
---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.159 on 293 degrees of freedom
Multiple R-squared: 0.6168, Adjusted R-squared: 0.6089
F-statistic: 78.59 on 6 and 293 DF, p-value: <2.2e-16
> plot(fit)


R2R^2R2是0.6左右不是很显著&＃xff0c;系数的$t$检验中常数项$x4$和$x5$都不显著。回归方程的检验显著。还有一个模型的检验图如下。

5.2模型分析

上图的分析方式参见文献³的第八章。

5.3方差分析表

在R语言里是可以输出方差分析表的&＃xff1a;

> anova(fit)
Analysis of Variance Table
下面准备剔除变量介绍一般的变量选择的方法。
Response: yDf Sum Sq Mean Sq F value Pr(>F)
x1 1 1508.74 1508.74 323.7227 <2.2e-16 ***
x2 1 363.08 363.08 77.9040 <2.2e-16 ***
x3 1 290.80 290.80 62.3947 5.704e-14 ***
x4 1 3.63 3.63 0.7792 0.37812
x5 1 0.40 0.40 0.0857 0.76995
x6 1 31.12 31.12 6.6766 0.01025 *
Residuals 293 1365.55 4.66
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

5、变量选择

在多元线性回归中变量选择的方法众多&＃xff0c;包含向前向后回归所有子集以及套索压缩系数法。本节以上面的数据模型为例演示这几种方式。

5.1 逐步回归

逐步回归的实现需要用到MASS包中的step或者stepAIC函数。
下面是step的具体参数

> step
function (object, scope, scale &＃61; 0, direction &＃61; c("both", "backward", "forward"), trace &＃61; 1, keep &＃61; NULL, steps &＃61; 1000, use.start &＃61; FALSE, k &＃61; 2, ...)


参数说明&＃xff1a;

• object&＃xff1a;模型的对象
• diretion&＃xff1a;选择方法&＃xff0c;both逐步回归&＃xff0c;backward向前&＃xff0c;forward向后。
  代码案例&＃xff1a;

> step(fit, direction &＃61; c("both")) # 逐步回归
Start: AIC&＃61;468.66
y ~ x1 &＃43; x2 &＃43; x3 &＃43; x4 &＃43; x5 &＃43; x6Df Sum of Sq RSS AIC
- x5 1 0.33 1365.9 466.73
- x4 1 2.28 1367.8 467.16
1365.5 468.66
- x6 1 31.12 1396.7 473.42
- x1 1 36.60 1402.2 474.59
- x3 1 274.85 1640.4 521.67
- x2 1 396.64 1762.2 543.16
Step: AIC&＃61;466.73
y ~ x1 &＃43; x2 &＃43; x3 &＃43; x4 &＃43; x6Df Sum of Sq RSS AIC
- x4 1 2.24 1368.1 465.23
1365.9 466.73
&＃43; x5 1 0.33 1365.5 468.66
- x6 1 31.18 1397.1 471.51
- x1 1 36.44 1402.3 472.63
- x3 1 281.26 1647.2 520.91
- x2 1 398.89 1764.8 541.60
Step: AIC&＃61;465.23
y ~ x1 &＃43; x2 &＃43; x3 &＃43; x6Df Sum of Sq RSS AIC
1368.1 465.23
&＃43; x4 1 2.24 1365.9 466.73
&＃43; x5 1 0.29 1367.8 467.16
- x6 1 32.57 1400.7 470.28
- x1 1 35.65 1403.8 470.94
- x3 1 282.58 1650.7 519.55
- x2 1 404.83 1773.0 540.99
Call:
lm(formula &＃61; y ~ x1 &＃43; x2 &＃43; x3 &＃43; x6, data &＃61; Data)
Coefficients:
(Intercept) x1 x2 x3 x6 -0.04082 1.12071 1.10693 1.00183 1.04557
> step(fit, direction &＃61; c("backward")) # 向前回归
Start: AIC&＃61;468.66
y ~ x1 &＃43; x2 &＃43; x3 &＃43; x4 &＃43; x5 &＃43; x6Df Sum of Sq RSS AIC
- x5 1 0.33 1365.9 466.73
- x4 1 2.28 1367.8 467.16
1365.5 468.66
- x6 1 31.12 1396.7 473.42
- x1 1 36.60 1402.2 474.59
- x3 1 274.85 1640.4 521.67
- x2 1 396.64 1762.2 543.16
Step: AIC&＃61;466.73
y ~ x1 &＃43; x2 &＃43; x3 &＃43; x4 &＃43; x6Df Sum of Sq RSS AIC
- x4 1 2.24 1368.1 465.23
1365.9 466.73
- x6 1 31.18 1397.1 471.51
- x1 1 36.44 1402.3 472.63
- x3 1 281.26 1647.2 520.91
- x2 1 398.89 1764.8 541.60
Step: AIC&＃61;465.23
y ~ x1 &＃43; x2 &＃43; x3 &＃43; x6Df Sum of Sq RSS AIC
1368.1 465.23
- x6 1 32.57 1400.7 470.28
- x1 1 35.65 1403.8 470.94
- x3 1 282.58 1650.7 519.55
- x2 1 404.83 1773.0 540.99
Call:
lm(formula &＃61; y ~ x1 &＃43; x2 &＃43; x3 &＃43; x6, data &＃61; Data)
Coefficients:
(Intercept) x1 x2 x3 x6 -0.04082 1.12071 1.10693 1.00183 1.04557
> step(fit, direction &＃61; c("forward")) # 向后回归
Start: AIC&＃61;468.66
y ~ x1 &＃43; x2 &＃43; x3 &＃43; x4 &＃43; x5 &＃43; x6
Call:
lm(formula &＃61; y ~ x1 &＃43; x2 &＃43; x3 &＃43; x4 &＃43; x5 &＃43; x6, data &＃61; Data)
Coefficients:
(Intercept) x1 x2 x3 x4 -0.04555 1.13727 1.09968 0.99557 -0.08748 x5 x6 0.03488 1.02470


其中逐步回归和向前回归选择出了重要的变量&＃xff0c;而且估计的参数效果比较好。向后回归没有选择出好的变量。当然不是说向后回归选择变量不好&＃xff0c;可能换个数据向前回归也会失效。下面介绍所有子集法&＃xff0c;从所有可能的自变量的集合中选出最优的模型&＃xff0c;缺点是计算量比较大。这里变量只有6个&＃xff0c;所有我们使用相关函数进行计算。

5.2所有子集法

实现所有子集所用的包是leaps&＃xff0c;打开下面的函数参数需要用到泛型函数的知识&＃xff0c;具体参见泛型函数的讲解和案例⁷。

plot(x, labels&＃61;obj$xnames, main&＃61;NULL, scale&＃61;c("bic", "Cp", "adjr2", "r2"),
col&＃61;gray(seq(0, 0.9, length &＃61; 10)),...)


• scale这个参数是调整模型的选择顺序。例子如下&＃xff1a;

leap <- regsubsets(y~.,data &＃61; Data)
class(leap)
par(mfrow &＃61; c(2,2))
plot(leap,scale&＃61;"bic")
plot(leap,scale&＃61;"Cp")
plot(leap,scale&＃61;"adjr2")
plot(leap,scale&＃61;"r2")
?plot.regsubsets


可以看出$bic、cp、adjr2$三个标准可以选择出好的模型。r2的第二个模型是最优的。

5.3套索法

这个方法不同于上述两种方式&＃xff0c;是文献⁵的作者首次在1996年提出的。在文献二第三章的第四节中有此方法的具体介绍。本节使用ncvreg包来对上述模型实现变量选择。
先介绍ncvreg和cv.ncvreg这两个函数。

> ncvreg
function (X, y, family &＃61; c("gaussian", "binomial", "poisson"), penalty &＃61; c("MCP", "SCAD", "lasso"), gamma &＃61; switch(penalty, SCAD &＃61; 3.7, 3), alpha &＃61; 1, lambda.min &＃61; ifelse(n > p, 0.001, 0.05), nlambda &＃61; 100, lambda, eps &＃61; 1e-04, max.iter &＃61; 10000, convex &＃61; TRUE, dfmax &＃61; p &＃43; 1, penalty.factor &＃61; rep(1, ncol(X)), warn &＃61; TRUE, returnX, ...) > cv.ncvregfunction (X, y, ..., cluster, nfolds &＃61; 10, seed, fold, returnY &＃61; FALSE, trace &＃61; FALSE)


• X因变量矩阵
• 响应变量向量
• family 广义线性模型的分布族这里用高斯即可。
• penalty惩罚函数这里用MCP和lasso即可。
  案例&＃xff1a;

> X <- as.matrix(Data[,1:6])
> Y <- Data[,7]
> fit1 <- cv.ncvreg(X,Y,penalty &＃61; "MCP", trace &＃61; T)
Starting CV fold #1
Starting CV fold #2
Starting CV fold #3
Starting CV fold #4
Starting CV fold #5
Starting CV fold #6
Starting CV fold #7
Starting CV fold #8
Starting CV fold #9
Starting CV fold #10
> coef(fit1)
(Intercept) x1 x2 x3 x4 x5
-0.04211208 1.12333734 1.10523022 1.00117800 -0.02642666 0.00000000 x6 1.04100244
> summary(fit1)
MCP-penalized linear regression with n&＃61;300, p&＃61;6
At minimum cross-validation error (lambda&＃61;0.0685):
-------------------------------------------------Nonzero coefficients: 5Cross-validation error (deviance): 4.81R-squared: 0.60Signal-to-noise ratio: 1.47Scale estimate (sigma): 2.192
> fit2 <- ncvreg(X,Y,penalty &＃61; "lasso")
> fit3 <- ncvreg(X, Y, penalty &＃61; "lasso", lambda &＃61; 0.2)
> coef(fit3)
(Intercept) x1 x2 x3 x4 x5
-0.06064812 1.03875142 0.90957439 0.80761377 0.00000000 0.00000000 x6 0.94863056


plot(fit1)
plot(fit2)


这篇知识讲了回归模型的基本一些方法除此之外你还可能会用到回归诊断和异常值的改进&＃xff0c;最后列举一写常用的函数。