如何用Sklearn线性模型获取多元线性回归P值及优化与可视化

1. 解决P值获取问题 #

sklearn的LinearRegression仅专注于预测任务，不提供P值等统计检验指标，这是你无法获取P值的核心原因。改用statsmodels库的OLS模型可解决此问题，同时能对齐Excel分析结果：

关键注意事项 #

• Excel回归工具默认包含截距项，statsmodels需要手动添加常数项才能与Excel结果对齐
• 确保分类特征（如Company_Location）已正确做哑变量编码，否则会导致结果偏差
• 数据形状报错通常是因为未正确处理输入维度（比如忘记添加常数项，或特征矩阵与因变量维度不匹配）

修正代码示例 #

import pandas as pd
import statsmodels.api as sm

def get_regression_stats():
    dfreg = pd.read_csv("dfreg.csv")
    # 定义自变量和因变量
    X = dfreg[['Large_size', 'Mid_Level', 'Senior_Level', 'Exec_Level', 'Company_Location']]
    y = dfreg['Salary_In_USD']
    # 添加截距项（匹配Excel默认设置）
    X = sm.add_constant(X)
    # 拟合OLS模型并输出完整统计结果
    model = sm.OLS(y, X).fit()
    print(model.summary())
    return model

运行后会输出包含P值（对应P>|t|列）、系数、R²等完整统计信息，结果可与Excel分析工具库对齐。

2. 递归移除高P值特征的实现 #

以下函数会自动迭代训练模型，每次移除P值最大且超过0.05的特征，直到所有特征P值≤0.05或无特征剩余：

import pandas as pd
import statsmodels.api as sm

def recursive_feature_filter():
    dfreg = pd.read_csv("dfreg.csv")
    remaining_features = ['Large_size', 'Mid_Level', 'Senior_Level', 'Exec_Level', 'Company_Location']
    y = dfreg['Salary_In_USD']
    
    while True:
        if not remaining_features:
            print("所有特征均被移除")
            break
            
        X = dfreg[remaining_features]
        X = sm.add_constant(X)
        model = sm.OLS(y, X).fit()
        
        # 提取除截距外的特征P值
        p_values = model.pvalues.drop('const')
        max_p = p_values.max()
        
        if max_p <= 0.05:
            print("所有特征P值均符合要求")
            print(model.summary())
            break
        
        # 移除P值最高的特征
        feature_to_remove = p_values.idxmax()
        print(f"移除特征：{feature_to_remove}，P值：{max_p:.4f}")
        remaining_features.remove(feature_to_remove)

3. 数据集最佳可视化方法 #

针对你的分类特征为主的数据集，推荐以下可视化方式：

（1）特征系数与置信区间图 #

直观展示每个特征对薪资的影响幅度及统计显著性：

import matplotlib.pyplot as plt
import seaborn as sns

model = get_regression_stats()
# 提取系数和95%置信区间
coef_df = pd.DataFrame({
    '特征': model.params.index[1:],  # 排除截距项
    '系数': model.params.values[1:],
    '置信下限': model.conf_int().iloc[1:, 0],
    '置信上限': model.conf_int().iloc[1:, 1]
})

plt.figure(figsize=(10,6))
sns.barplot(x='系数', y='特征', data=coef_df, palette='viridis')
plt.errorbar(x=coef_df['系数'], y=coef_df['特征'], 
             xerr=[coef_df['系数']-coef_df['置信下限'], coef_df['置信上限']-coef_df['系数']], 
             fmt='none', c='black')
plt.title('特征系数与95%置信区间')
plt.show()

（2）分类特征与薪资的箱线图 #

展示不同分类特征下的薪资分布（如不同职位级别、公司规模的薪资差异）：

plt.figure(figsize=(12,6))
sns.boxplot(x='Senior_Level', y='Salary_In_USD', data=dfreg)
plt.title('资深级别与薪资分布')
plt.show()

# 可替换为其他分类特征（如Mid_Level、Exec_Level、Large_size等）

（3）残差分析图 #

验证回归模型的假设（残差正态性、同方差性）：

model = get_regression_stats()
residuals = model.resid

plt.figure(figsize=(12,5))
# 残差分布直方图
plt.subplot(1,2,1)
sns.histplot(residuals, kde=True)
plt.title('残差分布直方图')

# 残差QQ图
plt.subplot(1,2,2)
sm.qqplot(residuals, line='45', ax=plt.gca())
plt.title('残差QQ图')
plt.show()

（4）实际值vs预测值散点图 #

评估模型的拟合效果：

model = get_regression_stats()
y_pred = model.predict(sm.add_constant(dfreg[model.params.index[1:]]))

plt.figure(figsize=(8,8))
plt.scatter(dfreg['Salary_In_USD'], y_pred, alpha=0.6)
plt.plot([dfreg['Salary_In_USD'].min(), dfreg['Salary_In_USD'].max()], 
         [dfreg['Salary_In_USD'].min(), dfreg['Salary_In_USD'].max()], 
         'r--')
plt.xlabel('实际薪资')
plt.ylabel('预测薪资')
plt.title('实际vs预测薪资')
plt.show()

内容的提问来源于stack exchange，提问作者DataHawk1119