Python案例如何用Statsmodels做回归分析

本文目录导读：

1. 基础线性回归案例
2. 实际数据分析案例
3. 高级回归分析案例
4. 交互作用和分类变量处理
5. 异常值检测和处理
6. 实用技巧和最佳实践

我来详细说明如何使用Statsmodels进行回归分析,包含多个实际案例。

基础线性回归案例

案例1：房价预测

import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
# 创建示例数据
np.random.seed(42)
n = 100
# 生成模拟房价数据
sqft = np.random.normal(1500, 300, n)  # 房屋面积
bedrooms = np.random.randint(1, 5, n)  # 卧室数量
age = np.random.randint(0, 30, n)  # 房龄
location_score = np.random.uniform(1, 10, n)  # 位置评分
# 生成目标变量（房价）
price = 50000 + 200*sqft + 20000*bedrooms - 1000*age + 10000*location_score + np.random.normal(0, 20000, n)
# 创建DataFrame
data = pd.DataFrame({
    'sqft': sqft,
    'bedrooms': bedrooms,
    'age': age,
    'location_score': location_score,
    'price': price
})
print("数据概览：")
print(data.head())
print(f"\n样本数：{len(data)}")
print(f"特征数：{len(data.columns) - 1}")
# 准备数据
X = data[['sqft', 'bedrooms', 'age', 'location_score']]
y = data['price']
# 添加常数项（截距）
X = sm.add_constant(X)
# 拆分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# 拟合模型
model = sm.OLS(y_train, X_train).fit()
# 输出结果
print("\n" + "="*50)
print("回归分析结果")
print("="*50)
print(model.summary())
# 预测
y_pred = model.predict(X_test)
# 评估模型
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error
print("\n" + "="*50)
print("模型评估指标")
print("="*50)
print(f"R² 分数：{r2_score(y_test, y_pred):.3f}")
print(f"RMSE：{np.sqrt(mean_squared_error(y_test, y_pred)):.2f}")
print(f"MAE：{mean_absolute_error(y_test, y_pred):.2f}")

案例2：多元线性回归诊断

# 回归诊断
print("\n" + "="*50)
print("回归诊断")
print("="*50)
# 1. 残差分析
residuals = model.resid
fitted_values = model.fittedvalues
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# 残差 vs 拟合值
axes[0, 0].scatter(fitted_values, residuals, alpha=0.6)
axes[0, 0].axhline(y=0, color='r', linestyle='--')
axes[0, 0].set_xlabel('Fitted Values')
axes[0, 0].set_ylabel('Residuals')
axes[0, 0].set_title('Residuals vs Fitted')
# Q-Q图
sm.qqplot(residuals, line='s', ax=axes[0, 1])
axes[0, 1].set_title('Q-Q Plot')
# 残差直方图
axes[1, 0].hist(residuals, bins=20, edgecolor='black', alpha=0.7)
axes[1, 0].set_xlabel('Residuals')
axes[1, 0].set_ylabel('Frequency')
axes[1, 0].set_title('Residuals Distribution')
# 残差 vs 预测值排序
axes[1, 1].scatter(range(len(residuals)), residuals, alpha=0.6)
axes[1, 1].axhline(y=0, color='r', linestyle='--')
axes[1, 1].set_xlabel('Observation Order')
axes[1, 1].set_ylabel('Residuals')
axes[1, 1].set_title('Residuals vs Order')
plt.tight_layout()
plt.show()
# 2. 多重共线性检测（VIF）
from statsmodels.stats.outliers_influence import variance_inflation_factor
vif_data = pd.DataFrame()
vif_data["Feature"] = X.columns[1:]  # 排除常数项
vif_data["VIF"] = [variance_inflation_factor(X.values, i+1) for i in range(len(X.columns)-1)]
print("\n方差膨胀因子（VIF）分析：")
print(vif_data)
# 3. 正态性检验
from scipy import stats
shapiro_stat, shapiro_p = stats.shapiro(residuals)
print(f"\nShapiro-Wilk正态检验统计量：{shapiro_stat:.3f}")
print(f"p值：{shapiro_p:.4f}")
print("残差正态性：", "满足" if shapiro_p > 0.05 else "不满足")

实际数据分析案例

案例3：使用真实数据集（波士顿房价）

# 使用scikit-learn的波士顿房价数据集
from sklearn.datasets import load_diabetes
# 加载糖尿病数据集
diabetes = load_diabetes()
df = pd.DataFrame(diabetes.data, columns=diabetes.feature_names)
df['target'] = diabetes.target
print("糖尿病数据集信息：")
print(f"样本数：{len(df)}")
print(f"特征数：{len(df.columns) - 1}")
print(f"\n特征名称：{list(df.columns[:-1])}")
# 相关性分析
correlation_matrix = df.corr()
plt.figure(figsize=(12, 8))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', center=0)'特征相关性热力图')
plt.tight_layout()
plt.show()
# 选择显著特征
# 选取与目标变量相关性较高的特征
target_corr = correlation_matrix['target'].drop('target')
significant_features = target_corr[abs(target_corr) > 0.2].index.tolist()
print(f"\n显著特征（|corr| > 0.2）：{significant_features}")
# 建立回归模型
X = df[significant_features]
y = df['target']
X = sm.add_constant(X)
# 拟合模型
model_2 = sm.OLS(y, X).fit()
print("\n" + "="*50)
print("糖尿病数据回归分析结果")
print("="*50)
print(model_2.summary())
# 特征重要性可视化
feature_importance = pd.DataFrame({
    'Feature': model_2.params.index,
    'Coefficient': model_2.params.values,
    'P-Value': model_2.pvalues.values
})
feature_importance = feature_importance.sort_values('Coefficient', key=abs, ascending=False)
plt.figure(figsize=(10, 6))
bars = plt.bar(feature_importance['Feature'], feature_importance['Coefficient'])
plt.xticks(rotation=45)
plt.ylabel('Coefficient Value')'特征重要性（系数大小）')
plt.grid(True, alpha=0.3)
# 标记显著性
for i, (bar, pval) in enumerate(zip(bars, feature_importance['P-Value'])):
    if pval < 0.05:
        bar.set_color('green')
    elif pval < 0.1:
        bar.set_color('orange')
    else:
        bar.set_color('red')
plt.tight_layout()
plt.show()

高级回归分析案例

案例4：多项式回归和交互作用

# 创建非线性关系数据
np.random.seed(42)
x = np.linspace(0, 10, 100)
y = 2 + 1.5*x - 0.3*x**2 + np.random.normal(0, 1, 100)
# 创建多项式特征
X_poly = pd.DataFrame({
    'x': x,
    'x_squared': x**2,
    'x_cubic': x**3
})
X_poly = sm.add_constant(X_poly)
# 拟合不同阶数的多项式
print("="*50)
print("多项式回归比较")
print("="*50)
# 线性模型
model_linear = sm.OLS(y, sm.add_constant(pd.DataFrame({'x': x}))).fit()
print("\n线性模型：")
print(f"R²：{model_linear.rsquared:.3f}")
# 二次模型
model_quad = sm.OLS(y, X_poly[['const', 'x', 'x_squared']]).fit()
print(f"二次模型 R²：{model_quad.rsquared:.3f}")
# 三次模型
model_cubic = sm.OLS(y, X_poly).fit()
print(f"三次模型 R²：{model_cubic.rsquared:.3f}")
# 可视化比较
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.scatter(x, y, alpha=0.6, label='Actual data')
x_sorted = np.sort(x)
plt.plot(x_sorted, model_linear.predict(sm.add_constant(pd.DataFrame({'x': x_sorted}))), 
         'r-', label='Linear', alpha=0.8)
plt.plot(x_sorted, model_quad.predict(X_poly.sort_values('x')[['const', 'x', 'x_squared']]), 
         'g-', label='Quadratic', alpha=0.8)
plt.plot(x_sorted, model_cubic.predict(X_poly.sort_values('x')), 
         'b-', label='Cubic', alpha=0.8)
plt.xlabel('x')
plt.ylabel('y')'多项式回归比较')
plt.legend()
plt.grid(True, alpha=0.3)
plt.subplot(1, 2, 2)
# 残差比较
residuals_dict = {
    'Linear': model_linear.resid,
    'Quadratic': model_quad.resid,
    'Cubic': model_cubic.resid
}
for name, resid in residuals_dict.items():
    plt.hist(resid, bins=15, alpha=0.5, label=name)
plt.xlabel('Residuals')
plt.ylabel('Frequency')'残差分布比较')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 模型选择指标
print("\n模型选择指标：")
models = [model_linear, model_quad, model_cubic]
for i, model in enumerate(models):
    print(f"\n{'线性' if i==0 else '二次' if i==1 else '三次'}模型：")
    print(f"  R²：{model.rsquared:.4f}")
    print(f"  Adjusted R²：{model.rsquared_adj:.4f}")
    print(f"  AIC：{model.aic:.2f}")
    print(f"  BIC：{model.bic:.2f}")

交互作用和分类变量处理

案例5：包含分类变量的回归

# 创建包含分类变量的数据
np.random.seed(42)
n = 200
# 生成数据
area = np.random.choice(['A', 'B', 'C'], n)
size = np.random.normal(100, 30, n)
condition = np.random.choice(['Good', 'Fair', 'Poor'], n, p=[0.3, 0.4, 0.3])
# 生成目标变量
price = 200000 + 500*size
price += np.where(area == 'A', 50000, np.where(area == 'B', 20000, 0))
price += np.where(condition == 'Good', 30000, np.where(condition == 'Fair', 10000, -10000))
price += np.random.normal(0, 15000, n)
# 创建DataFrame
df_cat = pd.DataFrame({
    'area': area,
    'size': size,
    'condition': condition,
    'price': price
})
print("包含分类变量的数据：")
print(df_cat.head())
# 处理分类变量（one-hot编码）
X_cat = pd.get_dummies(df_cat[['area', 'size', 'condition']], 
                       columns=['area', 'condition'], 
                       drop_first=True)
X_cat = sm.add_constant(X_cat)
y_cat = df_cat['price']
# 拟合模型
model_cat = sm.OLS(y_cat, X_cat).fit()
print("\n" + "="*50)
print("包含分类变量的回归结果")
print("="*50)
print(model_cat.summary())
# 分析分类变量的影响
print("\n分类变量影响分析：")
for var in ['area_B', 'area_C', 'condition_Fair', 'condition_Poor']:
    if var in model_cat.params.index:
        coef = model_cat.params[var]
        pval = model_cat.pvalues[var]
        print(f"{var}: 系数={coef:.2f}, p值={pval:.4f}")

异常值检测和处理

案例6：稳健回归

# 创建包含异常值的数据
np.random.seed(42)
x = np.linspace(0, 10, 50)
y = 2 + 1.5*x + np.random.normal(0, 1, 50)
# 添加异常值
y[45:] = y[45:] + 20
# 普通最小二乘法
X_rob = sm.add_constant(x)
ols_model = sm.OLS(y, X_rob).fit()
# 稳健回归（Huber）
huber_model = sm.RLM(y, X_rob, M=sm.robust.norms.HuberT()).fit()
# 比较结果
print("="*50)
print("OLS vs 稳健回归比较")
print("="*50)
print("\nOLS 系数：")
print(ols_model.params)
print("\n稳健回归（Huber）系数：")
print(huber_model.params)
# 可视化
plt.figure(figsize=(10, 6))
plt.scatter(x, y, alpha=0.6, label='Data')
plt.plot(x, ols_model.params[0] + ols_model.params[1]*x, 
         'r-', label='OLS', linewidth=2)
plt.plot(x, huber_model.params[0] + huber_model.params[1]*x, 
         'g--', label='Huber (Robust)', linewidth=2)
plt.xlabel('x')
plt.ylabel('y')'OLS vs 稳健回归比较')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
# 权重分析
print("\n异常值权重：")
print(f"最后一个点的权重：{huber_model.weights[-1]:.3f}")
print(f"正常点的权重：{huber_model.weights[0]:.3f}")

实用技巧和最佳实践

# 1. 特征选择 - 逐步回归
def stepwise_regression(X, y, threshold_in=0.01, threshold_out=0.05):
    """
    逐步向前回归
    """
    included = []
    while True:
        changed = False
        # 向前选入
        excluded = list(set(X.columns) - set(included))
        new_pval = pd.Series(index=excluded)
        for new_column in excluded:
            model = sm.OLS(y, sm.add_constant(X[included + [new_column]])).fit()
            new_pval[new_column] = model.pvalues[new_column]
        best_pval = new_pval.min()
        if best_pval < threshold_in:
            best_feature = new_pval.idxmin()
            included.append(best_feature)
            changed = True
            print(f"Forward Selection: Add {best_feature} with p-value {best_pval:.4f}")
        # 向后剔除
        while True:
            model = sm.OLS(y, sm.add_constant(X[included])).fit()
            # 排除常数项
            pvalues = model.pvalues.iloc[1:]
            worst_pval = pvalues.max()
            if worst_pval > threshold_out:
                worst_feature = pvalues.idxmax()
                included.remove(worst_feature)
                print(f"Backward Elimination: Remove {worst_feature} with p-value {worst_pval:.4f}")
            else:
                break
        if not changed:
            break
    return included
# 2. 交叉验证
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LinearRegression
# 使用sklearn的交叉验证
lr = LinearRegression()
cv_scores = cross_val_score(lr, X, y, cv=5, scoring='r2')
print(f"\n交叉验证 R² 分数：{cv_scores}")
print(f"平均 R²：{cv_scores.mean():.3f} (+/- {cv_scores.std() * 2:.3f})")
# 3. 模型解释
# 系数置信区间
print("\n系数置信区间：")
print(model.conf_int())

使用Statsmodels进行回归分析的关键步骤：

1. 数据准备：处理缺失值，编码分类变量，添加常数项
2. 模型建立：选择合适的回归方法（OLS, RLM, WLS等）
3. 模型诊断：检查残差、多重共线性、异常值
4. 模型改进：多项式项、交互作用、变量变换
5. 模型评估：R²、AIC/BIC、交叉验证
6. 结果解释：系数解读、显著性检验、置信区间

这些案例覆盖了Statsmodels回归分析的主要应用场景,你可以根据实际需求选择适合的方法。