Python 实战多元线性回归模型，附带原理+代码

作者 | 萝卜

来源 | 早起Python（ ID：zaoqi-python ）

「多元线性回归模型」非常常见，是大多数人入门机器学习的第一个案例，尽管如此，里面还是有许多值得学习和注意的地方。其中多元共线性这个问题将贯穿所有的机器学习模型，所以本文会「将原理知识穿插于代码段中」，争取以不一样的视角来叙述和讲解「如何更好的构建和优化多元线性回归模型」。主要将分为两个部分：

• 详细原理

• Python 实战

Python 实战

Python 多元线性回归的模型的实战案例有非常多，这里虽然选用的经典的房价预测，但贵在的流程简洁完整，其中用到的精度优化方法效果拔群，能提供比较好的参考价值。

数据探索

本文的数据集是经过清洗的美国某地区的房价数据集。

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
 
df = pd.read_csv('house_prices.csv')
df.info();df.head()

参数说明：

• neighborhood/area: 所属街区和面积

• bedrooms/bathrooms: 卧室和浴室

• style: 房屋样式

多元线性回归建模

现在我们直接构建多元线性回归模型：

from statsmodels.formula.api import ols
                      # 小写的 ols 函数才会自带截距项，OLS 则不会
         # 固定格式：因变量 ~ 自变量(+ 号连接)
lm = ols('price ~ area + bedrooms + bathrooms', data=df).fit()
lm.summary()

红框为我们关注的结果值，其中截距项 Intercept 的 P 值没有意义，可以不用管它。

模型优化

从上图可以看到，模型的精度较低，因为还有类别变量 neighborhood 和 style 没有完全利用。这里我们先查看一下类别变量的类别分布情况：

# 类别变量，又称为名义变量，nominal variables
nominal_vars = ['neighborhood', 'style']
 
for each in nominal_vars:
    print(each, ':')
    print(df[each].agg(['value_counts']).T)  # Pandas 骚操作
    # 直接 .value_counts().T 无法实现下面的效果
     ## 必须得 agg，而且里面的中括号 [] 也不能少
    print('='*35)

虚拟变量的设置

因为类别变量无法直接放入模型，这里需要转换一下，而多元线性回归模型中类别变量的转换最常用的方法之一便是将其转化成虚拟变量。原理其实非常简单，将无法直接用于建模的名义变量转换成可放入模型的虚拟变量的核心就短短八个字：「四散拆开，非此即彼」。下面用一个只有 4 行的微型数据集辅以说明。

从上表中，不难发现：

• 该名义变量有 n 类，就能拆分出 n 个虚拟变量。

• 巧妙的使用 0 和 1 来达到「用虚拟变量列代替原名义变量所在类别」。

接下来要做的就是将生成的虚拟变量们放入多元线性回归模型，但要注意的是：「转化后的虚拟变量们需要舍弃一个」，才能得到满秩矩阵。具体原因和有关线性代数的解释可以查看笔者打包好的论文，我们可以理解为，当该名义变量可划分为 n 类时，只需要 n-1 个虚拟变量就已足够获知所有信息了。该丢弃哪个，可根据实际情况来决定。

因此为原数据集的某名义变量添加虚拟变量的步骤为:

• 抽出希望转换的名义变量(一个或多个)

• pandas 的 get_dummies 函数

• 与原数据集横向拼接

注意虚拟变量设置成功后，需要与原来的数据集拼接，这样才能将其一起放进模型。

再次建模后，发现模型精度大大提升，但潜在的多元共线性问题也随之显现出来。

在解释模型中虚拟变量的系数之前，我们先消除模型中多元共线性的影响，因为在排除共线性后，模型中的各个自变量的系数又会改变，最终的多元线性回归模型的等式又会不一样。多重线性回归模型的主要假设之一是我们的预测变量(自变量)彼此不相关。我们希望预测变量(自变量)与反应变量(因变量)相关，而不是彼此之间具有相关性。方差膨胀因子 ( Variance Inflation Factor，以下简称 VIF )，是「指解释变量之间存在多重共线性时的方差与不存在多重共线性时的方差之比」。

上图公式可以看出在方差膨胀因子的检测中:

• 每个自变量都会有一个膨胀因子值  ，最后根据值的大小来选择是否删减。

• 既然   表示相关性，是谁跟谁的相关性呢? 是自变量中的某一变量与除它外剩余的自变量进行多元线性回归，取回归结果，即模型精度来作为这个变量与剩余自变量的相关性。听起来可能有点绕，这里举一下实例：用 “面积、卧室数量和浴室数量” 作为自变量来预测房价，在进行自变量的方差膨胀因子的检测时，面积、卧室数和浴室数轮流做单独的因变量，剩下的两个变量作为自变量，来看看这三个自变量中哪个变量对其余两个变量的解释性高。

•    越大，如已经到了 0.9，那分母就很小，  值就等于 10，即表示这个自变量已经同时解释了另外的某个或多个自变量，存在多元共线性，可以考虑删除一些自变量。

VIF 越大，显示共线性越严重。经验判断方法表明：当 0< VIF <10，不存在多重共线性；当10≤ VIF < 100，存在较强的多重共线性；当 VIF≥ 100，存在严重多重共线性。

方差膨胀因子的检测

我们自己来写一个方差膨胀因子的检测函数。

def vif(df, col_i):
    """
    df: 整份数据
    col_i：被检测的列名
    """
    cols = list(df.columns)
    cols.remove(col_i)
    cols_noti = cols
    formula = col_i + '~' + '+'.join(cols_noti)
    r2 = ols(formula, df).fit().rsquared
    return 1. / (1. - r2)

现在进行检测。

test_data = results[['area', 'bedrooms', 'bathrooms', 'A', 'B']]
for i in test_data.columns:
    print(i, '\t', vif(df=test_data, col_i=i))

发现 bedrooms 和 bathrooms 存在强相关性，可能这两个变量是解释同一个问题,方差膨胀因子较大的自变量通常是成对出现的。

果然，bedrooms 和 bathrooms 这两个变量的方差膨胀因子较高，这里删除自变量 bedrooms 再次进行建模。

lm = ols(formula='price ~ area + bathrooms + A + B', data=results).fit()
lm.summary()

模型精度稍降，但消除了多元共线性后能够使模型的泛化能力提升。再次进行多元共线性检测。

test_data = results[['area', 'bedrooms', 'A', 'B']]
for i in test_data.columns:
    print(i, '\t', vif(df=test_data, col_i=i))

那么多元共线性就只有通过方差膨胀因子才能看的出来吗？ 其实并不一定，通过结合散点图或相关稀疏矩阵和模型中自变量的系数也能看出端倪。下图是未处理多元共线性时的自变量系数。

可以很明显的看出，bathrooms 的参数很可能是有问题的，怎么可能 bathrooms 的数据量每增加一个，房屋总价还减少 1.373*10 的四次方美元呢？简单的画个散点图和热力图也应该知道房屋总价与 bathrooms 个数应该是成正比例关系的。

模型解释

多元线性回归模型的可解释性比较强，将模型参数打印出来即可求出因变量与自变量的关系。

所以最终的建模结果如下，且该模型的精度为0.916。

另外在等式结果中，截距项 Intercept 和 area, bedrooms 等变量的系数都还好理解；A，B 这两个虚拟变量可能相对困难些。其实根据原理部分的表格来看，如果房屋在 C 区，那等式中 A 和 B 这两个字母的值便是 0，所以这便引出了非常重要的一点：使用了虚拟变量的多元线性回归模型结果中，存在于模型内的虚拟变量都是跟被删除掉的那个虚拟变量进行比较。所以这个结果便表示在其他情况完全一样时(即除虚拟变量外的项) A 区的房屋比 C 区低 8707.18 美元，B 区则比 C 区贵 449896.73.7 美元。当然我们也可以画个箱线图来查看与检验，发现结果正如模型中 A 与 B 的系数那般显示。

本文以多元线性回归为基础和前提，在因变量房价与多个自变量的实际观测值建立了多元线性回归模型；分析并检验各个预测变量对因变量的综合线性影响的显著性，并尽可能的消除多重共线性的影响，筛选出因变量有显著线性影响的自变量，对基准模型进行优化，并对各自变量相对重要性进行评定，进而提升了回归模型的预测精度。

更多精彩推荐

☞大写“惨”？三次改变世界、却被无情出局的程序员大牛！
☞110+优秀作品、20+热门领域曝光， 高通如何玩转“AI 的夏天”！
☞加码 2000 亿还不够，阿里云再建三座超级数据中心加速新基建
☞我还没考试，算法就说我的物理一定挂科
☞中台架构详解（上） | 大咖说中台
☞Eth2 验证者快速启动器发布，还有什么惊喜是我们不知道的？
点分享点点赞点在看