上一篇文章中介绍了一元线性回归(R语言数据分析与挖掘(第四章):回归分析(1)——一元回归分析),然而,在实际操作中,多元性回归会更多见,因为一个响应变量会对应多个解释变量,一种现象常常是与多个因素相联系的,由多个自变量的最优组合共同来预测或估计因变量,比只用一个自变量进行预测或估计更有效,更符合实际。因此多元线性回归比一元线性回归的实用意义更大。
关于多元线性回归的模型,百度百科有介绍(https://baike.baidu.com/item/多元线性回归/10702248?fr=aladdin),我们这里是R语言,重点是介绍怎么使用R语言实现多元线性回归分析。关于多元线性回归的模型在第二章(R语言数据分析与挖掘(第二章):统计学基础(视频))是有介绍的,因为这些都是统计学的基础。所以这里就不介绍了。没有打好基础的同学,先停下来,不要急。
下面我们继续iris为例:
# install.packages("car") # 如果没有安装该包,自行安装> library(car)> scatterplotMatrix(iris,1:4)
对于变量Petal.Length而言,与其变量Sepal.Length和Petal.Width呈正相关关系,而与Sepal.Width呈负相关关系。这种相关性可以使用函数cor()量化出来。
> coor = cor(iris[,1:4])> coor Sepal.Length Sepal.Width Petal.Length Petal.WidthSepal.Length 1.0000000 -0.1175698 0.8717538 0.8179411Sepal.Width -0.1175698 1.0000000 -0.4284401 -0.3661259Petal.Length 0.8717538 -0.4284401 1.0000000 0.9628654Petal.Width 0.8179411 -0.3661259 0.9628654 1.0000000
根据前面的分析,将Petal.Length作为响应变量,其他三个变量作为解释变量进行多元线性回归是合理的选择。
> lm2 = lm(Petal.Length~.,data = iris[,1:4])> summary(lm2)
Call:lm(formula = Petal.Length ~ ., data = iris[, 1:4])
Residuals: Min 1Q Median 3Q Max -0.99333 -0.17656 -0.01004 0.18558 1.06909
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.26271 0.29741 -0.883 0.379 Sepal.Length 0.72914 0.05832 12.502 <2e-16 ***Sepal.Width -0.64601 0.06850 -9.431 <2e-16 ***Petal.Width 1.44679 0.06761 21.399 <2e-16 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.319 on 146 degrees of freedomMultiple R-squared: 0.968, Adjusted R-squared: 0.9674 F-statistic: 1473 on 3 and 146 DF, p-value: < 2.2e-16
lm()函数中,我们直接指定Petal.Length~.,"."是代表有多个变量。截距项(Intercept)的检验值为为0.379, 没有通过显著性检验,说明我们需要对模型进行修正。
> lm3 = lm(Petal.Length~. 0,data = iris[,1:4])> summary(lm3)
Call:lm(formula = Petal.Length ~ . 0, data = iris[, 1:4])
Residuals: Min 1Q Median 3Q Max -0.98066 -0.18621 -0.02385 0.18795 1.05950
Coefficients: Estimate Std. Error t value Pr(>|t|) Sepal.Length 0.69421 0.04284 16.20 <2e-16 ***Sepal.Width -0.67286 0.06134 -10.97 <2e-16 ***Petal.Width 1.46802 0.06315 23.25 <2e-16 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3187 on 147 degrees of freedomMultiple R-squared: 0.9942, Adjusted R-squared: 0.9941 F-statistic: 8426 on 3 and 147 DF, p-value: < 2.2e-16
我们将截距去掉可以看到,通过了显著性检验。
查看残差拟合图:
> plot(lm3)Hit <Return> to see next plot: Hit <Return> to see next plot: Hit <Return> to see next plot: Hit <Return> to see next plot:
模型lm3的回归诊断可用函数ncvTest()查看。
> ncvTest(lm3)Non-constant Variance Score Test Variance formula: ~ fitted.values Chisquare = 7.831356, Df = 1, p = 0.0051348
p值小于0.05,说明回归拟合值不存在同方差性。此外,还需要考虑变量的交互作用。如在变量Sepal.Length和Petal.Width的交互作用,则需要添加Sepal.Length*Petal.Width项。
> lm4 = lm(Petal.Length ~ Sepal.Length *Petal.Width Sepal.Width,data= iris[,1:4])> summary(lm4)
Call:lm(formula = Petal.Length ~ Sepal.Length * Petal.Width Sepal.Width, data = iris[, 1:4])
Residuals: Min 1Q Median 3Q Max -0.92738 -0.20315 -0.01323 0.20128 1.03535
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.14477 0.54103 -2.116 0.0361 * Sepal.Length 0.87530 0.09482 9.231 3.00e-16 ***Petal.Width 2.03589 0.31038 6.559 8.91e-10 ***Sepal.Width -0.61154 0.07013 -8.719 5.93e-15 ***Sepal.Length:Petal.Width -0.10424 0.05362 -1.944 0.0539 . ---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.316 on 145 degrees of freedomMultiple R-squared: 0.9688, Adjusted R-squared: 0.968 F-statistic: 1127 on 4 and 145 DF, p-value: < 2.2e-16
上述代码在模型Im3的基础上添加交互作用,得到新的回归模型Im4,模型的摘要显示:在0.1的显著性水平下,截距项、三个解释变量以及交互项的参数估计均通过了显著性检验,表明该交互项的添加有一定的合理性。
在交互项的选择方面,原则上需要将解释变量进行组合,建模并参考R-squared项进行选取,使得R-squared变大且参数估计能通过显著性检验的交互项就可以引入回归模型中,该方法适用于解释变量不多的情况,在实际操作中,往往需要根据行业知识来判断解释变量间的交互作用。
下面根据训练的模型进行后续分析,由上述探讨可知,回归模型Im3是有效的,可以利用其进行后续的数据探索,如进行预测。此处采用训练集和测试集均为iris 数据集,并将其应用一般需要将原始数据随机拆分为训练集(train data)和测试集即采用原始的数据集训练模型并进行预测工作,然而在实际操作中,训练模型(test data)。
> pred = predict(lm3,data = iris[,1:4])> length(pred)[1] 150> plot(pred)
上诉代码表示,将iris数据集的前4列带入lm3回归模型中,预测变量Petal.Length的值,输出结果是150的向量,绘制预测值的散点图。