【python】逐步回归（多元线性回归模型中的应用）

---

前言

• Matlab中逐步回归的实现可以使用 Matlab 的 stepwise 函数，本文主要讨论逐步回归如何在 python 中使用。
• 思路参考视频：
  https://www.bilibili.com/video/BV1kU4y1R7o2/?spm_id_from=333.1007.top_right_bar_window_history.content.click&vd_source=67471d3a1b4f517b7a7964093e62f7e6

一、逐步回归

• 逐步回归（Stepwise Regression）是一种选择统计模型的技术，用于找到最优模型，即通过添加或移除变量来选择合适的特征。
• 逐步回归主要有三种方法：前进法（Forward Selection）、后退法（Backward Elimination）和逐步回归法（Stepwise Regression）。下面是对这三种方法的简单介绍：

1. 前进法（Forward Selection）

概念：

• 前进法从一个空模型开始，即最初没有任何预测变量。
• 然后逐步添加预测变量，每次添加一个变量，使得新模型的评价指标（如AIC、BIC、R^2等）最优。
• 继续这个过程，直到添加任何更多的变量都不能显著提高模型的性能。

步骤：

1. 从空模型开始，不包含任何预测变量。
2. 评估每个未加入模型的变量，将使模型性能最优的变量加入模型。
3. 重复步骤2，直到添加任何变量都不能显著改善模型。

2. 后退法（Backward Elimination）

概念：

• 后退法从包含所有预测变量的模型开始。
• 然后逐步移除预测变量，每次移除一个变量，使得新模型的评价指标最优。
• 继续这个过程，直到移除任何更多的变量都不能显著提高模型的性能。

步骤：

1. 从包含所有可能的预测变量的全模型开始。
2. 评估每个变量的显著性，移除最不显著的变量（即对模型贡献最小的变量）。
3. 重复步骤2，直到移除任何变量都不能显著改善模型。

3. 逐步回归法（Stepwise Regression）

概念：

• 逐步回归法结合了前进法和后退法，既可以添加变量也可以移除变量。
• 每次步骤既可以是添加一个新变量，也可以是移除一个现有变量，以达到模型性能的最优。

步骤：

1. 从空模型开始或包含所有预测变量的模型开始（具体取决于实现方式）。
2. 在每一步中，评估所有可能的添加或移除变量的操作。
3. 选择对模型性能最优的操作（添加或移除一个变量）。
4. 重复步骤2和步骤3，直到添加或移除任何变量都不能显著改善模型。

二、示例

• 水泥凝固时放出的热量 $y$ 与水泥中 4 种化学成分 $x_1,x_2,x_3,x_4$ 有关，今测得一组数据如下,试用逐步回归确定一个线性模型，并找出影响水泥凝固时放出热量的必要因素
  
  根据此示例，本文选用后退法选择出影响水泥凝固时放出热量的必要因素。

三、代码实现----python

1. 输入数据

import pandas as pd
import numpy as np
# 数据
x1 = np.array([7, 1, 11, 11, 7, 11, 3, 1, 2, 21, 1, 11, 10])
x2 = np.array([26, 29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68])
x3 = np.array([6, 15, 8, 8, 6, 9, 17, 22, 18, 4, 23, 9, 8])
x4 = np.array([60, 52, 20, 47, 33, 22, 6, 44, 22, 26, 34, 12, 12])
y = np.array([78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7, 72.5, 93.1, 115.9, 83.8, 113.3, 109.4])

# 自变量矩阵
X = pd.DataFrame({'x1': x1, 'x2': x2, 'x3': x3, 'x4': x4})

2. 初始化

本文选用的后退法，所以被选择的因素初始化为包含所有因素，被排出的元素列表为空。

# 初始化未被选中的因素
excluded = list(initial_list)
# 初始化被选中的因素
included = list(set(X.columns) - set(excluded))

3. 评估每个变量的显著性

本文中使用 $P$ 值的大小评估每个变量的显著性。

步骤：

1. 拟合回归模型
2. 获取所有特征的P值
3. 找到最大的P值及其对应的特征。

# 拟合选中的因素
model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included]))).fit()

# 得出被选中的因素里P值最大的
p_values = model.pvalues

# 排除常数项（截距项）的P值
p_values = p_values.drop('const')

# 找到最大的P值及其对应的特征
max_p_value = p_values.max()
print("最大的P值为:",max_p_value)

4. 判断函数退出的标志

直到被选中的因素拟合后得到的最大 $P$ 值小于 $\alpha (0.05)$，意味着移除任何变量都不能显著改善模型，函数退出。

def stepwise_selection(X, y, 
                       initial_list = [], 
                       threshold = 0.05,
                       mark = True):
                       
    # 初始化未被选中的因素
    excluded = list(initial_list)
    # 初始化被选中的因素
    included = list(set(X.columns) - set(excluded))

    while mark:

        # 拟合选中的因素
        model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included]))).fit()

        # 得出被选中的因素里P值最大的
        p_values = model.pvalues

        # 排除常数项（截距项）的P值
        p_values = p_values.drop('const')

        # 找到最大的P值及其对应的特征
        max_p_value = p_values.max()
        print("最大的P值为:",max_p_value)
        if max_p_value < threshold:
            mark = False
            print("最终模型:")
            print(model.summary())
            # 获取回归系数
            b = model.params
        else:
            max_p_feature = p_values.idxmax()
            print("最大的P值对应的特征为:",max_p_feature)

            # 从被选中的因素中去除
            included.remove(max_p_feature)
            print("更新后的因素为:",included)

    return included, b

5. 逐步回归的完整代码

import pandas as pd
import numpy as np
import statsmodels.api as sm

# 数据
x1 = np.array([7, 1, 11, 11, 7, 11, 3, 1, 2, 21, 1, 11, 10])
x2 = np.array([26, 29, 56, 31, 52, 55, 71, 31, 54, 47, 40, 66, 68])
x3 = np.array([6, 15, 8, 8, 6, 9, 17, 22, 18, 4, 23, 9, 8])
x4 = np.array([60, 52, 20, 47, 33, 22, 6, 44, 22, 26, 34, 12, 12])
y = np.array([78.5, 74.3, 104.3, 87.6, 95.9, 109.2, 102.7, 72.5, 93.1, 115.9, 83.8, 113.3, 109.4])

# 自变量矩阵
X = pd.DataFrame({'x1': x1, 'x2': x2, 'x3': x3, 'x4': x4})


def stepwise_selection(X, y, 
                       initial_list=[], 
                       threshold=0.05,
                       mark = True):
        
    # 初始化未被选中的因素
    excluded = list(initial_list)
    # 初始化被选中的因素
    included = list(set(X.columns) - set(excluded))

    while mark:

        # 拟合选中的因素
        model = sm.OLS(y, sm.add_constant(pd.DataFrame(X[included]))).fit()

        # 得出被选中的因素里P值最大的
        p_values = model.pvalues

        # 排除常数项（截距项）的P值
        p_values = p_values.drop('const')

        # 找到最大的P值及其对应的特征
        max_p_value = p_values.max()
        print("最大的P值为:",max_p_value)
        if max_p_value < threshold:
            mark = False
            print("最终模型:")
            print(model.summary())
            # 获取回归系数
            b = model.params
        else:
            max_p_feature = p_values.idxmax()
            print("最大的P值对应的特征为:",max_p_feature)

            # 从被选中的因素中去除
            included.remove(max_p_feature)
            print("更新后的因素为:",included)

    return included, b

result, b = stepwise_selection(X, y)

print(result)
print(b)

运行结果：

评估变量显著性的过程：

最终的模型：

回归系数：

6. 生成三维图的代码

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

# 计算Z的值
Z = b['const'] + b['x1'] * x1 + b['x2'] * x2

# 创建一个新的3D绘图对象
fig = plt.figure()
'''
111 的具体含义是：
第一位 1:整个图形只有 1 行。
第二位 1:整个图形只有 1 列。
第三位 1:子图在这个 1x1 网格中的第 1 个位置。
'''
ax = fig.add_subplot(111, projection='3d')

# 绘制散点图
ax.scatter(x1, x2, Z, c='r', marker='o', label='Data Points')

# 创建网格以绘制曲面
x1_range = np.linspace(min(x1), max(x1), 100)
x2_range = np.linspace(min(x2), max(x2), 100)
x1_grid, x2_grid = np.meshgrid(x1_range, x2_range)

# 计算曲面上的Z值
Z_surface = b['const'] + b['x1'] * x1_grid + b['x2'] * x2_grid

# 绘制曲面图
'''
alpha

意义:alpha 参数用于设置曲面的透明度。
取值范围:alpha 的取值范围是 0 到 1 之间。
0 表示完全透明，即不可见。
1 表示完全不透明。
作用：通过调整 alpha 参数，你可以在同一视图中更好地叠加多个图形，使得它们不会完全遮挡对方。

cmap

意义:cmap 参数用于设置曲面的颜色映射(colormap)。
常见的颜色映射:viridis, plasma, inferno, magma, cividis, jet, rainbow, coolwarm, hot 等。
作用：颜色映射用于根据 Z 值来着色曲面，帮助更清晰地展示高度或强度的变化。
'''
ax.plot_surface(x1_grid, x2_grid, Z_surface, alpha=0.5, cmap='viridis')

# 设置标签
ax.set_xlabel('X1')
ax.set_ylabel('X2')
ax.set_zlabel('Z')

# 添加图例
ax.legend()

# 显示图形
plt.show()

运行结果：