以前的文章已经讲述了简单线性回归的概念和代码实现,如今来继续看看多元线性回归。所谓多元线性回归其实就是自变量的个数变多了,以前的简单线性回归方程能够表示为:$y=b_0 +bx$,那么在多元中则是$y=b_0+b_1x_1+b_2x_2+...+b_nx_n$。python
在正式使用多元线性回归以前,咱们先谈谈关于线性回归的几个前置条件,首先,在线性回归中有几个重要的假设以下所示:算法
也就是说咱们在使用线性回归时,要考虑到实际状况是否吻合上述的几种假设,不然可能会致使拟合的模型不许确。api
除此以外,以前数据预处理的时候也提到过一个叫作虚拟变量的概念,如今对其进行详细的讲解。来看下面的数据集(一样因为篇幅问题只提供部分):app
R&D Spend | Administration | Marketing Spend | State | Profit |
---|---|---|---|---|
165349.2 | 136897.8 | 471784.1 | New York | 192261.83 |
162597.7 | 151377.59 | 443898.53 | California | 191792.06 |
153441.51 | 101145.55 | 407934.54 | Florida | 191050.39 |
144372.41 | 118671.85 | 383199.62 | New York | 182901.99 |
这组数据集反映的是公司的研发投入,行政支出,市场支出以及所在地点对于公司的收益的影响。其中所在地点这个变量是个分类变量,没有数值的概念,所以没法对其进行排序或者带入方程。数据预处理中,咱们使用了虚拟变量对其进行从新编码,那么对于这组数据咱们也能够作一样的操做,也就是说State这列数据能够拆分红三列:dom
New York | California | Florida |
---|---|---|
1 | 0 | 0 |
0 | 1 | 0 |
0 | 0 | 1 |
1 | 0 | 0 |
假设R&D Spend对应的自变量为$x_1$,Administration对应的是$x_2$,Marketing Spend对应的是$x_3$,虚拟编码对应的是$D_1$,$D_2$,$D_3$。可是对于虚拟编码,任意一个编码都能用另外两个来表示,好比是不是Florida,那么若是New York和California若是存在1,则确定不是Florida,不然就是Florida。那么能够定义其多元线性回归方程为:函数
$$ y = b_0 + b_1*x_1 + b_2*x_2 + b_3*x_3 + b_4*D_1 + b_5*D_2 $$工具
这里顺便一提,咱们这里给予虚拟变量的值是0和1,那么若是是(100,0)或者(-1,1)能够吗?答案是确定的,好比说(100,0),那咱们在方程中将前面的参数D缩小一百倍,达到的效果实际上是同样的。测试
接下来想一想,若是把咱们删掉的$b_6*D_3$加上那么这个方程对不对呢?其实这里就是虚拟变量的一个陷阱。回过头看看上面所说的关于回归问题的几个假设,第五条,无多重共线性这个性质咱们是否知足了?任意一个变量实际上均可以用其余的变量来表示,好比$D_3=1-D_1-D_2$,那么就不符合无多重共线性这个要求了。所以,在使用虚拟变量的时候必定要注意实际上使用的虚拟变量的数量应该是始终都要忽略掉一个的。优化
接下来咱们再谈谈如何一步一步地去构建一个多元线性回归模型。在实际应用中,每每会遇到对于一个因变量y,有不少的自变量x1,x2等等,但这些自变量不是全部的都是对这个y的预测颇有帮助因素,咱们须要从其中剔除掉无用的元素来获得最合适的模型,那么此时就有一个问题,如何来选择这些自变量呢?这里有五种方法来创建模型:编码
其中的2,3,4三种方式,也是咱们最经常使用的,叫作Step Regression也就是逐步回归,它们的算法是相似的,只是应用的顺序有点不一样。接下来一种种的介绍这几种方法。
所谓All-in,就是把全部咱们认为的自变量全扔进去,通常有这几种状况使用:
All-in这个方法很简单,但通常都是一些特殊状况下或者一些外力干涉时才会使用,这种方法不做推荐。
反向淘汰这个算法的精髓在于对于每个这个模型的自变量来讲,他对咱们的模型的预测结果实际上是有影响力的,用统计学上的一个概念叫作P-value来形容它的影响力。在这里咱们先定义它的这个影响力是否显著的一个门槛也就是一个significance level(SL),先定义为0.05。接着第二步使用全部的自变量来拟合出一个模型。第三步,对于这个模型当中的每个自变量都来计算它的P值(P-value),来显示它对咱们模型有多大的影响力,而后咱们取这个最高的P值,假设这个P>SL,就继续往第四步,不然就算法结束。那么第四步,最高的P值对应的那个自变量咱们就要将它从咱们的模型中去除。第五步,去除了一个自变量后,在用剩下的自变量从新对模型进行拟合,所以这里就是一个第三步到第五步的一个循环,直到全部剩下的P值都比SL要小,这样就说明模型已经拟合好了。步骤详情以下:
顺向选择和反向淘汰的算法思想很接近,只是顺序有了个颠倒。这里第一步仍是先定一个显著性的门槛SL=0.05.第二步,咱们在这边对每一个自变量$x_n$都进行简单线性回归拟合,分别获得它们的P值,而后取得它们中最低的。第三步对于这个最低的P值,咱们的结论就是这个自变量它对咱们将要拟合的模型的影响是最大的,因此说,咱们会保留这个自变量。接下来第四步咱们再看剩下的自变量当中加上哪个会给咱们带来最小的P值,假如说新的P值比咱们以前定义的SL小,那就从新回到第三步,也就是第三步又加入了一个新的变量而后在接下来剩下的变量当中,从新找最大的P值,而后继续加到模型当中,以此类推,直到剩下来的尚未被加入模型当中的变量它们的P值所有大于SL,也就是说剩下的变量对于模型有可能产生的影响不够显著,那对这些就不予采纳,这个时候模型就拟合好了。步骤详情以下:
所谓双向淘汰,其实就是对以前的两种算法的结合。在第一步中,咱们须要选择两个显著性门槛:一个旧的变量是否应该被剔除和一个新的尚未被采纳的变量是否应当进入咱们的模型。在第二步中,咱们要进行顺向选择,来决定是否采纳一个新的自变量。第三步要进行反向淘汰,也就是咱们可能要剔除旧的变量,而后在第二第三步之间进行循环,因为已经定义了两个门槛,但出现新的出不去,旧的进不来时,就说明模型已经拟合好了。步骤详情以下:
最后一种,信息量比较。所谓信息量,就是对于一个多元线性回归模型的一个评价方式。好比说给它进行打分。那么就有不少种打分方式。好比最多见的一种,叫作Akaike criterion.对于全部可能的模型,咱们对它们进行逐一的打分,对于多元线性回归,若是有N个自变量,那么就有$2^N-1$个不一样的模型,对这些模型打分后选择分数最高的模型。那这里就会有一个问题,若是N很大的时候,模型的数量就会很是庞大。因此说这个方法虽然直觉上很好理解,但自变量数量很大时就不适合使用这种方法。
接下来就要进行代码编写了,首先仍是老样子,先进行数据预处理,将数据导入并对分类数据进行虚拟编码,而后将数据集划分红训练集和测试集。上面有说起过虚拟编码的陷阱,实际上咱们使用的包中已经避开了这个陷阱,但为了强调这个问题,这里也对其进行处理一下:
import pandas as pd from sklearn.model_selection import train_test_split from sklearn.preprocessing import LabelEncoder, OneHotEncoder from sklearn.linear_model import LinearRegression import statsmodels.formula.api as sm import numpy as np data_path = '../data/50_Startups.csv' dataset = pd.read_csv(data_path) X = dataset.iloc[:, :-1].values y = dataset.iloc[:, 4].values labelencoder_X = LabelEncoder() X[:, 3] = labelencoder_X.fit_transform(X[:, 3]) onehotencoder = OneHotEncoder(categorical_features=[3]) X = onehotencoder.fit_transform(X).toarray() # Avoiding the Dummy Variable Trap X = X[:, 1:] X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)
预处理结束后,就来进行拟合线性回归模型,这里第一步拟合线性回归模型的方式和以前的简单线性回归是同样的。
# Fitting Multiple Linear Regression to the Training set regressor = LinearRegression() regressor.fit(X_train, y_train) # Predicting the Test set results y_pred = regressor.predict(X_test)
目前咱们拟合的模型使用的方法是上述提到的All-in方法,但这个方法并非最好的,这里咱们使用Backward Elimination来对回归器进行优化,使用的工具是statsmodels.formula.api。再而后咱们还要作一件事情,在多元线性回归时,不少时候模型里面都会有一个常数$b_0$,至关于$b_0*1$。但在用到的标准库的函数里面是不包含这个常数项的,所以咱们须要在包含自变量的矩阵中加上一列,这一列所有都是1,$b_0$就是它的系数。代码以下:
X_train = np.append(arr=np.ones((40, 1)), values=X_train, axis=1)
接下来要真正开始反向淘汰了,首先建立一个矩阵X_opt,包含了最佳的自变量选择,第一步其实是全部的自变量,以后的不断循环后会渐渐淘汰其中的自变量。
X_opt = X_train[:, [0, 1, 2, 3, 4, 5]]
这里为何用X_train[:, [0, 1, 2, 3, 4, 5]]而不是直接X_train呢?这是由于后面会进行不断的淘汰,会不断删除其中的列,后续的代码就能看出来了。
建立完X_opt后,接下来就要用它来拟合新的回归器regressor_OLS。这里用的sm.OLS方法的参数解释一下:endog对应的是因变量,这里就是y_train,exog指的是自变量,所以这里就是X_opt。
regressor_OLS = sm.OLS(endog=y_train, exog=X_opt).fit() regressor_OLS.summary()
而后使用summary()方法,这个方法能给咱们提供不少回归器的信息,执行后获得结果以下:
这里能看到全部自变量对用的P值,其中最大的是$x_2$,就是这个公司是否在加利福利亚洲,而后咱们剔除$x_2$,再继续拟合,获得summary,根据结果不断剔除自变量直到最终的P值都小于咱们定义的SL=0.05.代码以下:
X_opt = X_train[:, [0, 1, 3, 4, 5]] regressor_OLS = sm.OLS(endog=y_train, exog=X_opt).fit() regressor_OLS.summary() X_opt = X_train[:, [0, 3, 4, 5]] regressor_OLS = sm.OLS(endog=y_train, exog=X_opt).fit() regressor_OLS.summary() X_opt = X_train[:, [0, 3, 5]] regressor_OLS = sm.OLS(endog=y_train, exog=X_opt).fit() regressor_OLS.summary() X_opt = X_train[:, [0, 3]] regressor_OLS = sm.OLS(endog=y_train, exog=X_opt).fit() regressor_OLS.summary()
最终获得的结果是:
那么根据这个结果能够获得的结论就是一家公司的收益主要跟公司的研发投入有关。固然其实其中也有其余的自变量对应了很低的P值,若是本身一步步运行这段代码会发现行政的投入对应的P值也只有0.07不算很高,如何更好的判断线性回归模型的优劣后面还会有其余的方法来判断。