未记录事实核验 · 状态：维护中

记录仅说明核验时间与维护状态，不保证适用于所有环境；使用前请留意正文中的版本、适用条件与风险提示。

warning 适用边界与维护说明

本文适合用于模拟住房数据上的残差图、影响点筛查及模型处理策略学习，不宜直接用作生产模型验收或假设成立的认证。当前普通KS检验使用同批残差估计分布参数，其p值不能按参数预先已知的普通KS检验解释；正态性与同方差性检验不显著不能证明假设成立，也不能以多数票或任一项不显著作确认。检查清单仍将拟合值与绝对残差的相关性用于线性关系判定、以DW经验区间确认独立性，两处VIF辅助回归未加入常数项，异常点计数还存在重叠相加问题，因此各项结论及自动建议需另行审查。住房记录没有有意义的时间顺序，自相关部分仅供接口演示，前十阶最小p值不能作为一次预先指定检验的结果。Goldfeld-Quandt示例当前使用等大小、相同设计项的子组，不应未经调整推广至不同残差自由度的分组。Cook距离、杠杆值和VIF阈值仅为筛查线索，不应据此自动删除观测；不同样本或响应尺度的训练指标不能直接比较策略优劣，预测效果需在一致目标尺度及独立验证数据上评估。本次仅作全文结构审阅，未执行代码或开展外部事实核验。

🩺 回归诊断

abstract 📋 本节概览

🎯 学习目标：全面掌握回归分析中的诊断方法，能够识别和处理模型假设违反、异常值、影响点等问题
📚 前置知识：回归分析|回归分析 | 线性回归|线性回归 | 相关分析|相关分析 | 假设检验流程|假设检验流程
⏱️ 预计用时：40分钟
🛠️ 工具准备：Python、NumPy、SciPy、Pandas、Matplotlib、Statsmodels、Seaborn、Scikit-learn

note 📖 **核心定义**

回归诊断（Regression Diagnostics）：用于评估回归模型的合理性、检验模型假设、识别异常值和影响点、评估模型预测能力的一系列统计方法和图形技术。

目的：确保回归模型的可靠性和有效性，为模型改进和决策提供依据。

📐 回归模型的基本假设

线性回归的核心假设

note 🔍 **回归分析的五大假设**

条件均值与线性设定：模型对参数线性，条件均值由所选自变量及变换项表示，并要求E[误差|自变量]=0；这不等于每个原始自变量与因变量都呈直线关系。
误差相关结构：常规独立样本推断需要相应的独立性条件；时间序列或聚类数据需要适配其误差相关结构，不能仅凭残差图或DW经验区间确认独立性。
同方差性：经典同方差标准误要求误差的条件方差恒定；异方差不必然使OLS系数失去有效性，但通常需要调整标准误或估计方法。
正态性：条件误差正态性主要支持经典有限样本精确推断，并非拟合线性回归或所有大样本推断的统一必要条件；拟合残差只是误差的诊断代理。
无完全多重共线性：设计矩阵需满列秩才能唯一估计系数；高但非完全共线性主要影响估计精度与稳定性，不是禁止拟合的统一条件。

违反假设的后果

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error, r2_score
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan, het_white
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.stattools import durbin_watson
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体和绘图样式
plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS']
plt.rcParams['axes.unicode_minus'] = False
sns.set_style("whitegrid")

def demonstrate_assumption_violations():
    """
    演示违反回归假设的影响
    """
    print("=== 回归假设违反的影响演示 ===\n")
    
    np.random.seed(42)
    n = 100
    
 # 生成基础数据
    x = np.random.uniform(0, 10, n)
    
 # 1. 正常情况（满足所有假设）
    y_normal = 2 + 3 * x + np.random.normal(0, 2, n)
    
 # 2. 非线性关系
    y_nonlinear = 2 + 3 * x + 0.5 * x**2 + np.random.normal(0, 2, n)
    
 # 3. 异方差（方差随x增大）
    y_hetero = 2 + 3 * x + np.random.normal(0, 0.5 + 0.5*x, n)
    
 # 4. 非正态残差（使用t分布）
    y_nonnormal = 2 + 3 * x + stats.t.rvs(df=3, size=n) * 2
    
    cases = {
        '满足假设': y_normal,
        '非线性关系': y_nonlinear, 
        '异方差': y_hetero,
        '非正态残差': y_nonnormal
    }
    
    print("不同假设违反情况下的回归结果比较:")
    print(f"{'情况':<12} {'R²':<8} {'回归系数':<10} {'标准误':<10} {'p值':<12}")
    print("-" * 60)
    
    for case_name, y in cases.items():
 # 拟合线性回归
        X = sm.add_constant(x)
        model = sm.OLS(y, X).fit()
        
        coef = model.params[1]  # x的系数
        se = model.bse[1]       # 标准误
        p_value = model.pvalues[1]  # p值
        r2 = model.rsquared
        
        print(f"{case_name:<12} {r2:<8.3f} {coef:<10.3f} {se:<10.3f} {p_value:<12.6f}")
    
    print("\n📊 影响分析:")
    print("• 非线性关系: 降低拟合优度，系数估计偏差")
    print("• 异方差: 标准误估计不准确，假设检验失效")
    print("• 非正态残差: 置信区间和假设检验可能不可靠")
    
    return x, cases

# 执行假设违反演示
x_demo, y_cases = demonstrate_assumption_violations()

🔍 残差分析

残差的类型和计算

def comprehensive_residual_analysis():
    """
    全面的残差分析
    """
    print("\n=== 全面残差分析 ===")
    
 # 生成示例数据：房价预测
    np.random.seed(123)
    n = 150
    
 # 自变量：面积、楼层、房龄
    area = np.random.normal(100, 30, n)
    floor = np.random.randint(1, 31, n)
    age = np.random.exponential(5, n)
    
 # 因变量：房价（万元），含一些异常值
    price = (2.5 * area + 0.3 * floor - 0.8 * age + 50 + 
             np.random.normal(0, 20, n))
    
 # 添加几个异常值
    outlier_indices = [10, 25, 87, 134]
    for idx in outlier_indices:
        price[idx] += np.random.choice([-80, 80])  # 价格异常
    
 # 创建数据框
    df = pd.DataFrame({
        'area': area,
        'floor': floor, 
        'age': age,
        'price': price
    })
    
    print("房价预测模型数据概况:")
    print(df.describe().round(2))
    
 # 拟合多元线性回归
    X = df[['area', 'floor', 'age']]
    y = df['price']
    
 # 添加常数项
    X_with_const = sm.add_constant(X)
    model = sm.OLS(y, X_with_const).fit()
    
    print(f"\n模型拟合结果:")
    print(model.summary())
    
 # 计算各种残差
    fitted_values = model.fittedvalues
    residuals = model.resid  # 普通残差
    standardized_residuals = model.resid_pearson  # 标准化残差
    studentized_residuals = model.get_influence().resid_studentized_internal  # 学生化残差
    
 # 计算其他诊断统计量
    influence = model.get_influence()
    leverage = influence.hat_matrix_diag  # 杠杆值
    cooks_distance = influence.cooks_distance[0]  # Cook距离
    
    print(f"\n残差统计摘要:")
    print(f"残差均值: {np.mean(residuals):.6f}")
    print(f"残差标准差: {np.std(residuals):.3f}")
    print(f"最大标准化残差: {np.max(np.abs(standardized_residuals)):.3f}")
    print(f"最大杠杆值: {np.max(leverage):.3f}")
    print(f"最大Cook距离: {np.max(cooks_distance):.3f}")
    
 # 识别异常点
    print(f"\n异常点识别:")
    
 # 1. 高残差点（|标准化残差| > 2）
    high_residual_idx = np.where(np.abs(standardized_residuals) > 2)[0]
    print(f"高残差点 (|标准化残差|>2): {len(high_residual_idx)}个")
    if len(high_residual_idx) > 0:
        print(f"   索引: {high_residual_idx}")
    
 # 2. 高杠杆点（杠杆值 > 2(p+1)/n）
    p = X.shape[1]  # 自变量个数
    leverage_threshold = 2 * (p + 1) / n
    high_leverage_idx = np.where(leverage > leverage_threshold)[0]
    print(f"高杠杆点 (杠杆值>{leverage_threshold:.3f}): {len(high_leverage_idx)}个")
    if len(high_leverage_idx) > 0:
        print(f"   索引: {high_leverage_idx}")
    
 # 3. 高影响点（Cook距离 > 4/n）
    cook_threshold = 4 / n
    high_influence_idx = np.where(cooks_distance > cook_threshold)[0]
    print(f"高影响点 (Cook距离>{cook_threshold:.3f}): {len(high_influence_idx)}个")
    if len(high_influence_idx) > 0:
        print(f"   索引: {high_influence_idx}")
    
    return df, model, {
        'fitted': fitted_values,
        'residuals': residuals,
        'std_residuals': standardized_residuals,
        'stud_residuals': studentized_residuals,
        'leverage': leverage,
        'cooks_d': cooks_distance
    }

# 执行残差分析
df_housing, housing_model, residual_data = comprehensive_residual_analysis()

残差图和诊断图

def create_diagnostic_plots(df, model, residual_data):
    """
    创建回归诊断图
    """
    fig, axes = plt.subplots(3, 3, figsize=(20, 18))
    fig.suptitle('回归诊断图集', fontsize=16, fontweight='bold')
    
    fitted = residual_data['fitted']
    residuals = residual_data['residuals']
    std_residuals = residual_data['std_residuals']
    leverage = residual_data['leverage']
    cooks_d = residual_data['cooks_d']
    
 # 1. 残差vs拟合值图（检查线性关系和同方差性）
    axes[0, 0].scatter(fitted, residuals, alpha=0.6)
    axes[0, 0].axhline(y=0, color='red', linestyle='--', alpha=0.8)
    axes[0, 0].set_xlabel('拟合值')
    axes[0, 0].set_ylabel('残差')
    axes[0, 0].set_title('残差 vs 拟合值')
    
 # 添加局部回归线
    from scipy.interpolate import make_interp_spline
    if len(fitted) > 10:
        sorted_idx = np.argsort(fitted)
        fitted_sorted = fitted[sorted_idx]
        residuals_sorted = residuals[sorted_idx]
        
        try:
 # 使用移动平均来绘制趋势线
            window = max(10, len(fitted)//10)
            trend = pd.Series(residuals_sorted).rolling(window, center=True).mean()
            axes[0, 0].plot(fitted_sorted, trend, color='blue', linewidth=2, alpha=0.8)
        except:
            pass
    
    axes[0, 0].grid(True, alpha=0.3)
    
 # 2. Q-Q图（检查正态性）
    stats.probplot(residuals, dist="norm", plot=axes[0, 1])
    axes[0, 1].set_title('残差正态性Q-Q图')
    axes[0, 1].grid(True, alpha=0.3)
    
 # 3. 标准化残差vs拟合值（检查异方差）
    axes[0, 2].scatter(fitted, np.abs(std_residuals), alpha=0.6)
    axes[0, 2].set_xlabel('拟合值')
    axes[0, 2].set_ylabel('|标准化残差|')
    axes[0, 2].set_title('标准化残差绝对值 vs 拟合值')
    axes[0, 2].grid(True, alpha=0.3)
    
 # 4. 杠杆值vs标准化残差（识别影响点）
    axes[1, 0].scatter(leverage, std_residuals, alpha=0.6)
    axes[1, 0].set_xlabel('杠杆值')
    axes[1, 0].set_ylabel('标准化残差')
    axes[1, 0].set_title('杠杆值 vs 标准化残差')
    
 # 添加参考线
    axes[1, 0].axhline(y=2, color='red', linestyle='--', alpha=0.5)
    axes[1, 0].axhline(y=-2, color='red', linestyle='--', alpha=0.5)
    p = df.shape[1] - 1  # 自变量个数
    leverage_threshold = 2 * (p + 1) / len(df)
    axes[1, 0].axvline(x=leverage_threshold, color='red', linestyle='--', alpha=0.5)
    axes[1, 0].grid(True, alpha=0.3)
    
 # 5. Cook距离图
    indices = range(len(cooks_d))
    axes[1, 1].stem(indices, cooks_d, basefmt=" ")
    cook_threshold = 4 / len(df)
    axes[1, 1].axhline(y=cook_threshold, color='red', linestyle='--', alpha=0.8, 
                      label=f'阈值 = {cook_threshold:.3f}')
    axes[1, 1].set_xlabel('观测点索引')
    axes[1, 1].set_ylabel('Cook距离')
    axes[1, 1].set_title('Cook距离图')
    axes[1, 1].legend()
    axes[1, 1].grid(True, alpha=0.3)
    
 # 6. 残差直方图
    axes[1, 2].hist(residuals, bins=20, density=True, alpha=0.7, color='skyblue', edgecolor='black')
    
 # 叠加正态分布曲线
    mu, sigma = stats.norm.fit(residuals)
    x_norm = np.linspace(residuals.min(), residuals.max(), 100)
    axes[1, 2].plot(x_norm, stats.norm.pdf(x_norm, mu, sigma), 'r-', linewidth=2, label='正态分布拟合')
    axes[1, 2].set_xlabel('残差值')
    axes[1, 2].set_ylabel('密度')
    axes[1, 2].set_title('残差分布直方图')
    axes[1, 2].legend()
    axes[1, 2].grid(True, alpha=0.3)
    
 # 7. 按自变量的残差图（以面积为例）
    if 'area' in df.columns:
        axes[2, 0].scatter(df['area'], residuals, alpha=0.6)
        axes[2, 0].axhline(y=0, color='red', linestyle='--', alpha=0.8)
        axes[2, 0].set_xlabel('面积')
        axes[2, 0].set_ylabel('残差')
        axes[2, 0].set_title('残差 vs 面积')
        axes[2, 0].grid(True, alpha=0.3)
    
 # 8. 拟合值vs观测值
    y_true = model.model.endog
    axes[2, 1].scatter(fitted, y_true, alpha=0.6)
    
 # 添加y=x参考线
    min_val = min(fitted.min(), y_true.min())
    max_val = max(fitted.max(), y_true.max())
    axes[2, 1].plot([min_val, max_val], [min_val, max_val], 'r--', alpha=0.8, label='y=x')
    axes[2, 1].set_xlabel('拟合值')
    axes[2, 1].set_ylabel('观测值')
    axes[2, 1].set_title('拟合值 vs 观测值')
    axes[2, 1].legend()
    axes[2, 1].grid(True, alpha=0.3)
    
 # 9. 影响点综合图（Cook距离 vs 杠杆值）
    scatter = axes[2, 2].scatter(leverage, cooks_d, c=np.abs(std_residuals), 
                                cmap='viridis', alpha=0.6, s=50)
    axes[2, 2].set_xlabel('杠杆值')
    axes[2, 2].set_ylabel('Cook距离')
    axes[2, 2].set_title('影响点诊断图')
    
 # 添加阈值线
    axes[2, 2].axhline(y=4/len(df), color='red', linestyle='--', alpha=0.5, label='Cook阈值')
    axes[2, 2].axvline(x=2*(p+1)/len(df), color='red', linestyle='--', alpha=0.5, label='杠杆阈值')
    axes[2, 2].legend()
    axes[2, 2].grid(True, alpha=0.3)
    
 # 添加颜色条
    plt.colorbar(scatter, ax=axes[2, 2], label='|标准化残差|')
    
    plt.tight_layout()
    plt.show()
    
 # 诊断解释
    print("\n=== 诊断图解读指南 ===")
    print("📊 各图表的诊断要点:")
    print("1. 残差vs拟合值: 无明显模式表示未见明显设定偏离；漏斗形提示可能存在异方差，均不能单独证明假设成立")
    print("2. Q-Q图: 点接近直线表示与正态分布形状较一致；明显偏离提示进一步检查，不能单凭该图确认正态性")
    print("3. |标准化残差|vs拟合值: 随机散布→同方差；递增趋势→异方差")
    print("4. 杠杆值vs标准化残差: 识别高杠杆点和异常值")
    print("5. Cook距离: >4/n的点可列为潜在影响点，需结合敏感性分析和业务背景判断")
    print("6. 残差直方图: 检查残差分布的正态性")

# 创建诊断图
create_diagnostic_plots(df_housing, housing_model, residual_data)

🔗 多重共线性诊断

VIF方差膨胀因子

def multicollinearity_diagnosis():
    """
    多重共线性诊断
    """
    print("\n=== 多重共线性诊断 ===")
    
 # 创建包含共线性的数据
    np.random.seed(456)
    n = 200
    
 # 生成高度相关的自变量
    x1 = np.random.normal(0, 1, n)
    x2 = 0.8 * x1 + np.random.normal(0, 0.6, n)  # 与x1高度相关
    x3 = np.random.normal(0, 1, n)
    x4 = 2 * x1 + 3 * x2 + np.random.normal(0, 0.1, n)  # 几乎完全由x1和x2决定
    x5 = np.random.normal(0, 1, n)
    
 # 因变量
    y = 2 + 1.5*x1 + 2*x2 + 0.5*x3 - 1*x4 + 0.8*x5 + np.random.normal(0, 1, n)
    
 # 创建数据框
    df_mc = pd.DataFrame({
        'x1': x1, 'x2': x2, 'x3': x3, 'x4': x4, 'x5': x5, 'y': y
    })
    
    print("变量间相关系数矩阵:")
    correlation_matrix = df_mc.corr()
    print(correlation_matrix.round(3))
    
 # 计算VIF
    X = df_mc[['x1', 'x2', 'x3', 'x4', 'x5']]
    
    print(f"\n方差膨胀因子(VIF):")
    print(f"{'变量':<8} {'VIF':<10} {'诊断结果':<20}")
    print("-" * 40)
    
    vif_data = []
    for i in range(X.shape[1]):
        vif = variance_inflation_factor(X.values, i)
        vif_data.append((X.columns[i], vif))
        
 # VIF诊断标准
        if vif < 5:
            diagnosis = "无共线性问题"
        elif vif < 10:
            diagnosis = "中等共线性"
        else:
            diagnosis = "严重共线性"
            
        print(f"{X.columns[i]:<8} {vif:<10.2f} {diagnosis:<20}")
    
 # 特征值诊断
    print(f"\n相关矩阵特征值诊断:")
    eigenvalues = np.linalg.eigvals(X.corr())
    condition_indices = np.sqrt(eigenvalues.max() / eigenvalues)
    
    print(f"{'特征值':<10} {'条件指数':<10} {'诊断':<15}")
    print("-" * 35)
    
    for i, (eigenval, ci) in enumerate(zip(eigenvalues, condition_indices)):
        if ci < 10:
            diagnosis = "无问题"
        elif ci < 30:
            diagnosis = "中度共线性"
        else:
            diagnosis = "严重共线性"
            
        print(f"{eigenval:<10.4f} {ci:<10.2f} {diagnosis:<15}")
    
 # 展示共线性影响
    print(f"\n共线性对回归结果的影响:")
    
 # 完整模型
    X_full = sm.add_constant(X)
    model_full = sm.OLS(y, X_full).fit()
    
 # 移除高VIF变量后的模型
    X_reduced = X[['x1', 'x3', 'x5']]  # 移除x2和x4
    X_reduced_const = sm.add_constant(X_reduced)
    model_reduced = sm.OLS(y, X_reduced_const).fit()
    
    print(f"\n模型比较:")
    print(f"{'模型':<15} {'R²':<8} {'调整R²':<10} {'AIC':<10}")
    print("-" * 45)
    print(f"{'完整模型':<15} {model_full.rsquared:<8.3f} {model_full.rsquared_adj:<10.3f} {model_full.aic:<10.1f}")
    print(f"{'简化模型':<15} {model_reduced.rsquared:<8.3f} {model_reduced.rsquared_adj:<10.3f} {model_reduced.aic:<10.1f}")
    
 # 可视化相关性矩阵
    plt.figure(figsize=(10, 8))
    mask = np.triu(np.ones_like(correlation_matrix, dtype=bool))
    sns.heatmap(correlation_matrix, mask=mask, annot=True, cmap='coolwarm', 
                center=0, square=True, linewidths=0.5, cbar_kws={"shrink": .8})
    plt.title('变量相关性热图')
    plt.tight_layout()
    plt.show()
    
    return df_mc, vif_data, model_full, model_reduced

# 执行多重共线性诊断
df_multicollin, vif_results, full_model, reduced_model = multicollinearity_diagnosis()

📊 假设检验与统计测试

正态性检验

def normality_tests(residuals):
    """
    残差正态性检验
    """
    print("\n=== 残差正态性检验 ===")
    
 # 1. Shapiro-Wilk检验
    shapiro_stat, shapiro_p = stats.shapiro(residuals)
    print(f"Shapiro-Wilk检验:")
    print(f"  统计量: {shapiro_stat:.4f}")
    print(f"  p值: {shapiro_p:.6f}")
    print(f"  结论: {'拒绝正态性' if shapiro_p < 0.05 else '接受正态性'} (α=0.05)")
    
 # 2. Anderson-Darling检验
    ad_stat, ad_critical, ad_significance = stats.anderson(residuals, dist='norm')
    print(f"\nAnderson-Darling检验:")
    print(f"  统计量: {ad_stat:.4f}")
    
    for i, (sig, crit) in enumerate(zip(ad_significance, ad_critical)):
        result = "拒绝正态性" if ad_stat > crit else "未拒绝正态性原假设"
        print(f"  显著性水平 {sig}%: 临界值={crit:.3f}, {result}")
    
 # 3. Kolmogorov-Smirnov检验
    ks_stat, ks_p = stats.kstest(residuals, 'norm', args=(np.mean(residuals), np.std(residuals)))
    print(f"\nKolmogorov-Smirnov检验:")
    print(f"  统计量: {ks_stat:.4f}")
    print(f"  p值: {ks_p:.6f}")
    print(f"  结论: {'拒绝正态性' if ks_p < 0.05 else '接受正态性'} (α=0.05)")
    
 # 4. Jarque-Bera检验
    jb_stat, jb_p = stats.jarque_bera(residuals)
    print(f"\nJarque-Bera检验:")
    print(f"  统计量: {jb_stat:.4f}")
    print(f"  p值: {jb_p:.6f}")
    print(f"  结论: {'拒绝正态性' if jb_p < 0.05 else '接受正态性'} (α=0.05)")
    
 # 正态性综合评价
    tests_results = [shapiro_p > 0.05, ks_p > 0.05, jb_p > 0.05]
    normal_count = sum(tests_results)
    
    print(f"\n📊 正态性综合评价:")
    print(f"通过正态性检验的测试: {normal_count}/3")
    
    if normal_count >= 2:
        print("✅ 残差基本满足正态性假设")
    elif normal_count == 1:
        print("⚠️ 残差正态性存在疑问，建议进一步检查")
    else:
        print("❌ 残差明显违反正态性假设，需要模型调整")

# 对住房模型进行正态性检验
normality_tests(residual_data['residuals'])

同方差性检验

def heteroscedasticity_tests(model, df):
    """
    异方差检验
    """
    print("\n=== 异方差检验 ===")
    
 # 1. Breusch-Pagan检验
    bp_stat, bp_p, bp_f_stat, bp_f_p = het_breuschpagan(model.resid, model.model.exog)
    print(f"Breusch-Pagan检验:")
    print(f"  LM统计量: {bp_stat:.4f}")
    print(f"  p值: {bp_p:.6f}")
    print(f"  结论: {'存在异方差' if bp_p < 0.05 else '同方差'} (α=0.05)")
    
 # 2. White检验
    white_stat, white_p, white_f_stat, white_f_p = het_white(model.resid, model.model.exog)
    print(f"\nWhite检验:")
    print(f"  LM统计量: {white_stat:.4f}")
    print(f"  p值: {white_p:.6f}")
    print(f"  结论: {'存在异方差' if white_p < 0.05 else '同方差'} (α=0.05)")
    
 # 3. 戈德菲尔德-匡特检验（针对特定变量）
    if 'area' in df.columns:
 # 按面积排序后分组检验
        df_sorted = df.sort_values('area')
        n = len(df_sorted)
        n1 = n // 3
        n2 = n - n1
        
 # 前1/3和后1/3的方差比较
        group1 = df_sorted.iloc[:n1]
        group2 = df_sorted.iloc[-n1:]
        
 # 分别拟合回归
        X1 = sm.add_constant(group1[['area', 'floor', 'age']])
        X2 = sm.add_constant(group2[['area', 'floor', 'age']])
        
        model1 = sm.OLS(group1['price'], X1).fit()
        model2 = sm.OLS(group2['price'], X2).fit()
        
 # 计算F统计量（方差比）
        var1 = model1.ssr / model1.df_resid
        var2 = model2.ssr / model2.df_resid
        f_stat = var2 / var1 if var2 > var1 else var1 / var2
        
 # 计算p值
        df1 = model1.df_resid
        df2 = model2.df_resid
        p_value = 2 * (1 - stats.f.cdf(f_stat, max(df1, df2), min(df1, df2)))
        
        print(f"\nGoldfeld-Quandt检验 (基于面积):")
        print(f"  F统计量: {f_stat:.4f}")
        print(f"  p值: {p_value:.6f}")
        print(f"  结论: {'存在异方差' if p_value < 0.05 else '同方差'} (α=0.05)")
    
 # 异方差性综合评价
    tests_results = [bp_p < 0.05, white_p < 0.05]
    hetero_count = sum(tests_results)
    
    print(f"\n📊 异方差性综合评价:")
    print(f"检测出异方差的测试: {hetero_count}/2")
    
    if hetero_count == 0:
        print("✅ 数据满足同方差性假设")
    elif hetero_count == 1:
        print("⚠️ 异方差性存在疑问，建议进一步检查")
    else:
        print("❌ 数据明显存在异方差，需要处理")
        print("💡 建议处理方法:")
        print("   • 使用异方差稳健标准误（HC标准误）")
        print("   • 对因变量进行变换（如对数变换）")
        print("   • 使用加权最小二乘法(WLS)")

# 对住房模型进行异方差检验
heteroscedasticity_tests(housing_model, df_housing)

自相关检验

def autocorrelation_tests(model):
    """
    自相关检验（主要用于时间序列数据）
    """
    print("\n=== 自相关检验 ===")
    
 # Durbin-Watson检验
    dw_stat = durbin_watson(model.resid)
    print(f"Durbin-Watson检验:")
    print(f"  DW统计量: {dw_stat:.4f}")
    
 # DW统计量解释
    if 1.5 < dw_stat < 2.5:
        dw_conclusion = "DW处于经验参考区间，不能据此确认无显著自相关"
    elif dw_stat < 1.5:
        dw_conclusion = "可能存在正自相关"
    else:
        dw_conclusion = "可能存在负自相关"
    
    print(f"  结论: {dw_conclusion}")
    print(f"  解释: DW≈2表示无自相关，DW<2表示正相关，DW>2表示负相关")
    
 # Ljung-Box检验
    from statsmodels.stats.diagnostic import acorr_ljungbox
    
    lb_results = acorr_ljungbox(model.resid, lags=10, return_df=True)
    print(f"\nLjung-Box检验 (前10阶滞后):")
    print(f"  最低p值: {lb_results['lb_pvalue'].min():.6f}")
    
    if lb_results['lb_pvalue'].min() < 0.05:
        print("  结论: 存在显著自相关")
    else:
        print("  结论: 无显著自相关")

# 对住房模型进行自相关检验
autocorrelation_tests(housing_model)

🛠️ 模型改进策略

异常值处理

def outlier_treatment_strategies():
    """
    异常值处理策略演示
    """
    print("\n=== 异常值处理策略 ===")
    
 # 使用之前的住房数据
    df_original = df_housing.copy()
    
 # 识别异常值
    influence = housing_model.get_influence()
    cooks_d = influence.cooks_distance[0]
    leverage = influence.hat_matrix_diag
    std_residuals = housing_model.resid_pearson
    
 # 定义异常值标准
    cook_threshold = 4 / len(df_original)
    leverage_threshold = 2 * (3 + 1) / len(df_original)  # 3个自变量
    residual_threshold = 2
    
 # 综合异常值识别
    outliers_cook = cooks_d > cook_threshold
    outliers_leverage = leverage > leverage_threshold
    outliers_residual = np.abs(std_residuals) > residual_threshold
    
 # 组合标准：满足任意两个条件
    outliers_combined = (outliers_cook.astype(int) + 
                        outliers_leverage.astype(int) + 
                        outliers_residual.astype(int)) >= 2
    
    outlier_indices = np.where(outliers_combined)[0]
    
    print(f"识别出的异常值数量: {len(outlier_indices)}")
    print(f"异常值索引: {outlier_indices}")
    
 # 策略1: 移除异常值
    df_removed = df_original.drop(df_original.index[outlier_indices])
    X_removed = sm.add_constant(df_removed[['area', 'floor', 'age']])
    model_removed = sm.OLS(df_removed['price'], X_removed).fit()
    
 # 策略2: Winsorizing（缩尾处理）
    from scipy.stats import mstats
    df_winsor = df_original.copy()
    
 # 对连续变量进行winsorizing
    for col in ['area', 'price']:
        df_winsor[col] = mstats.winsorize(df_winsor[col], limits=[0.02, 0.02])
    
    X_winsor = sm.add_constant(df_winsor[['area', 'floor', 'age']])
    model_winsor = sm.OLS(df_winsor['price'], X_winsor).fit()
    
 # 策略3: 稳健回归
    from sklearn.linear_model import HuberRegressor
    
    X_robust = df_original[['area', 'floor', 'age']]
    y_robust = df_original['price']
    
    huber_reg = HuberRegressor(epsilon=1.35, max_iter=100)
    huber_reg.fit(X_robust, y_robust)
    
 # 比较结果
    print(f"\n不同处理策略的效果比较:")
    print(f"{'策略':<15} {'R²':<8} {'RMSE':<10} {'样本量':<8}")
    print("-" * 45)
    
 # 原始模型
    y_pred_orig = housing_model.fittedvalues
    rmse_orig = np.sqrt(mean_squared_error(df_original['price'], y_pred_orig))
    print(f"{'原始模型':<15} {housing_model.rsquared:<8.3f} {rmse_orig:<10.2f} {len(df_original):<8}")
    
 # 移除异常值模型
    y_pred_removed = model_removed.fittedvalues
    rmse_removed = np.sqrt(mean_squared_error(df_removed['price'], y_pred_removed))
    print(f"{'移除异常值':<15} {model_removed.rsquared:<8.3f} {rmse_removed:<10.2f} {len(df_removed):<8}")
    
 # Winsorizing模型
    y_pred_winsor = model_winsor.fittedvalues
    rmse_winsor = np.sqrt(mean_squared_error(df_winsor['price'], y_pred_winsor))
    print(f"{'缩尾处理':<15} {model_winsor.rsquared:<8.3f} {rmse_winsor:<10.2f} {len(df_winsor):<8}")
    
 # 稳健回归
    y_pred_robust = huber_reg.predict(X_robust)
    r2_robust = r2_score(y_robust, y_pred_robust)
    rmse_robust = np.sqrt(mean_squared_error(y_robust, y_pred_robust))
    print(f"{'稳健回归':<15} {r2_robust:<8.3f} {rmse_robust:<10.2f} {len(df_original):<8}")
    
    print(f"\n💡 异常值处理建议:")
    print("• 移除异常值: 提高拟合度但可能丢失重要信息")
    print("• 缩尾处理: 保留所有样本，减少极端值影响")
    print("• 稳健回归: 自动降低异常值权重，推荐用于探索性分析")
    
    return {
        'original': housing_model,
        'removed': model_removed,
        'winsorized': model_winsor,
        'robust': huber_reg
    }

# 执行异常值处理策略比较
treatment_results = outlier_treatment_strategies()

模型变换和改进

def model_transformation_strategies():
    """
    模型变换策略
    """
    print("\n=== 模型变换策略 ===")
    
    df_transform = df_housing.copy()
    
 # 1. 对数变换（处理异方差和偏斜分布）
    print("1. 对数变换策略:")
    
 # 检查对数变换输入：price + 1和area必须为有限正数
    if (not np.isfinite(df_transform[['price', 'area']].to_numpy()).all()
            or (df_transform['price'] <= -1).any()
            or (df_transform['area'] <= 0).any()):
        raise ValueError('对数变换要求price和area为有限值，且price > -1、area > 0')
    df_transform['log_price'] = np.log(df_transform['price'] + 1)
    df_transform['log_area'] = np.log(df_transform['area'])
    
 # 拟合对数模型
    X_log = sm.add_constant(df_transform[['log_area', 'floor', 'age']])
    model_log = sm.OLS(df_transform['log_price'], X_log).fit()
    
    print(f"   对数模型 R²: {model_log.rsquared:.3f}")
    
 # 检查对数变换后的异方差性
    bp_stat_log, bp_p_log, _, _ = het_breuschpagan(model_log.resid, model_log.model.exog)
    print(f"   对数模型异方差检验 p值: {bp_p_log:.4f}")
    
 # 2. 多项式变换
    print(f"\n2. 多项式变换策略:")
    
 # 添加二次项
    df_transform['area_squared'] = df_transform['area'] ** 2
    df_transform['age_squared'] = df_transform['age'] ** 2
    
    X_poly = sm.add_constant(df_transform[['area', 'area_squared', 'floor', 'age', 'age_squared']])
    model_poly = sm.OLS(df_transform['price'], X_poly).fit()
    
    print(f"   多项式模型 R²: {model_poly.rsquared:.3f}")
    print(f"   调整R²: {model_poly.rsquared_adj:.3f}")
    
 # 3. 交互项
    print(f"\n3. 交互效应策略:")
    
 # 添加交互项
    df_transform['area_floor'] = df_transform['area'] * df_transform['floor']
    
    X_interaction = sm.add_constant(df_transform[['area', 'floor', 'age', 'area_floor']])
    model_interaction = sm.OLS(df_transform['price'], X_interaction).fit()
    
    print(f"   交互模型 R²: {model_interaction.rsquared:.3f}")
    
 # 4. Box-Cox变换
    print(f"\n4. Box-Cox变换:")
    
    from scipy.stats import boxcox
    
 # 确保所有价格为正
    price_positive = df_transform['price'] - df_transform['price'].min() + 1
    price_boxcox, lambda_optimal = boxcox(price_positive)
    
    print(f"   最优λ参数: {lambda_optimal:.3f}")
    
    df_transform['price_boxcox'] = price_boxcox
    X_boxcox = sm.add_constant(df_transform[['area', 'floor', 'age']])
    model_boxcox = sm.OLS(df_transform['price_boxcox'], X_boxcox).fit()
    
    print(f"   Box-Cox模型 R²: {model_boxcox.rsquared:.3f}")
    
 # 检查Box-Cox变换后的正态性
    shapiro_stat_bc, shapiro_p_bc = stats.shapiro(model_boxcox.resid)
    print(f"   Box-Cox残差正态性检验 p值: {shapiro_p_bc:.6f}")
    
 # 模型比较总结
    print(f"\n📊 模型变换效果对比:")
    print(f"{'模型类型':<15} {'R²':<8} {'调整R²':<8} {'AIC':<10} {'BIC':<10}")
    print("-" * 55)
    
    models_comparison = [
        ('原始模型', housing_model),
        ('对数模型', model_log),
        ('多项式模型', model_poly),
        ('交互模型', model_interaction),
        ('Box-Cox模型', model_boxcox)
    ]
    
    for name, model in models_comparison:
        print(f"{name:<15} {model.rsquared:<8.3f} {model.rsquared_adj:<8.3f} {model.aic:<10.1f} {model.bic:<10.1f}")
    
    print(f"\n💡 模型选择建议:")
    print("• 对相同观测、相同响应尺度且似然定义可比的模型，AIC/BIC通常越小越好；本表不能直接跨响应变换排名")
    print("• 调整R²考虑了模型复杂度")
    print("• 需要平衡拟合度和模型简洁性")
    print("• 检查残差假设是否得到改善")
    
    return models_comparison

# 执行模型变换策略
transformation_results = model_transformation_strategies()

📋 回归诊断检查清单

def regression_diagnostic_checklist(model, df, X_cols, y_col):
    """
    回归诊断完整检查清单
    """
    print("=" * 60)
    print("📋 回归诊断完整检查清单")
    print("=" * 60)
    
    results = {}
    
 # 1. 模型基本信息
    print("\n1️⃣ 模型基本信息")
    print("-" * 30)
    print(f"样本量: {len(df)}")
    print(f"自变量数量: {len(X_cols)}")
    print(f"R²: {model.rsquared:.4f}")
    print(f"调整R²: {model.rsquared_adj:.4f}")
    print(f"F统计量: {model.fvalue:.4f} (p={model.f_pvalue:.6f})")
    
 # 2. 线性关系检查
    print("\n2️⃣ 线性关系检查")
    print("-" * 30)
    
 # 检查残差vs拟合值的模式
    fitted = model.fittedvalues
    residuals = model.resid
    
 # 计算残差的非线性模式
    correlation_resid_fitted = stats.pearsonr(fitted, np.abs(residuals))[1]
    
    if correlation_resid_fitted > 0.05:
        linearity_result = "✅ 线性关系假设基本满足"
    else:
        linearity_result = "⚠️ 可能存在非线性关系"
    
    print(f"残差-拟合值相关性检验: {linearity_result}")
    results['linearity'] = correlation_resid_fitted > 0.05
    
 # 3. 独立性检查
    print("\n3️⃣ 独立性检查")
    print("-" * 30)
    
    dw_stat = durbin_watson(residuals)
    independence_result = "✅ 基本满足独立性" if 1.5 < dw_stat < 2.5 else "⚠️ 可能存在自相关"
    print(f"Durbin-Watson统计量: {dw_stat:.4f} - {independence_result}")
    results['independence'] = 1.5 < dw_stat < 2.5
    
 # 4. 同方差性检查
    print("\n4️⃣ 同方差性检查")
    print("-" * 30)
    
    bp_stat, bp_p, _, _ = het_breuschpagan(residuals, model.model.exog)
    white_stat, white_p, _, _ = het_white(residuals, model.model.exog)
    
    homoscedasticity_bp = bp_p > 0.05
    homoscedasticity_white = white_p > 0.05
    homoscedasticity_overall = homoscedasticity_bp and homoscedasticity_white
    
    print(f"Breusch-Pagan检验: p={bp_p:.6f} {'✅' if homoscedasticity_bp else '❌'}")
    print(f"White检验: p={white_p:.6f} {'✅' if homoscedasticity_white else '❌'}")
    print(f"综合结论: {'✅ 满足同方差性' if homoscedasticity_overall else '❌ 存在异方差'}")
    results['homoscedasticity'] = homoscedasticity_overall
    
 # 5. 正态性检查
    print("\n5️⃣ 正态性检查")
    print("-" * 30)
    
    shapiro_stat, shapiro_p = stats.shapiro(residuals)
    jb_stat, jb_p = stats.jarque_bera(residuals)
    
    normality_shapiro = shapiro_p > 0.05
    normality_jb = jb_p > 0.05
    normality_overall = normality_shapiro or normality_jb  # 至少一个通过
    
    print(f"Shapiro-Wilk检验: p={shapiro_p:.6f} {'✅' if normality_shapiro else '❌'}")
    print(f"Jarque-Bera检验: p={jb_p:.6f} {'✅' if normality_jb else '❌'}")
    print(f"综合结论: {'✅ 基本满足正态性' if normality_overall else '❌ 违反正态性'}")
    results['normality'] = normality_overall
    
 # 6. 多重共线性检查
    print("\n6️⃣ 多重共线性检查")
    print("-" * 30)
    
    X = df[X_cols]
    vif_values = []
    vif_names = []
    
    for i in range(X.shape[1]):
        vif = variance_inflation_factor(X.values, i)
        vif_values.append(vif)
        vif_names.append(X.columns[i])
        
    max_vif = max(vif_values)
    multicollinearity_result = max_vif < 5
    
    print(f"最大VIF值: {max_vif:.2f}")
    for name, vif in zip(vif_names, vif_values):
        status = "✅" if vif < 5 else "⚠️" if vif < 10 else "❌"
        print(f"  {name}: {vif:.2f} {status}")
    
    print(f"综合结论: {'✅ 无严重共线性' if multicollinearity_result else '❌ 存在多重共线性'}")
    results['multicollinearity'] = multicollinearity_result
    
 # 7. 异常值和影响点检查
    print("\n7️⃣ 异常值和影响点检查")
    print("-" * 30)
    
    influence = model.get_influence()
    cooks_d = influence.cooks_distance[0]
    leverage = influence.hat_matrix_diag
    std_residuals = model.resid_pearson
    
 # 计算阈值
    cook_threshold = 4 / len(df)
    leverage_threshold = 2 * (len(X_cols) + 1) / len(df)
    
 # 统计异常点
    high_cook = np.sum(cooks_d > cook_threshold)
    high_leverage = np.sum(leverage > leverage_threshold)
    high_residual = np.sum(np.abs(std_residuals) > 2)
    
    print(f"高Cook距离点 (>{cook_threshold:.3f}): {high_cook}个")
    print(f"高杠杆点 (>{leverage_threshold:.3f}): {high_leverage}个")
    print(f"高残差点 (|标准化残差|>2): {high_residual}个")
    
    outliers_acceptable = (high_cook + high_leverage + high_residual) < len(df) * 0.05  # 少于5%
    print(f"综合结论: {'✅ 异常点数量可接受' if outliers_acceptable else '⚠️ 异常点较多，需要关注'}")
    results['outliers'] = outliers_acceptable
    
 # 8. 总体评估
    print("\n8️⃣ 总体模型评估")
    print("-" * 30)
    
    passed_assumptions = sum([
        results['linearity'],
        results['independence'], 
        results['homoscedasticity'],
        results['normality'],
        results['multicollinearity'],
        results['outliers']
    ])
    
    total_assumptions = 6
    pass_rate = passed_assumptions / total_assumptions
    
    print(f"通过的假设检验: {passed_assumptions}/{total_assumptions} ({pass_rate*100:.1f}%)")
    
    if pass_rate >= 0.8:
        overall_conclusion = "示例规则下诊断通过率较高；不代表模型可靠，仍需结合业务与独立验证评估"
        color = "green"
    elif pass_rate >= 0.6:
        overall_conclusion = "示例规则下部分诊断触发警示；是否可用需结合用途和具体问题判断"
        color = "yellow"
    else:
        overall_conclusion = "示例规则下较多诊断触发警示；不能仅凭通过率认定必须重新建模"
        color = "red"
    
    print(f"\n🎯 最终结论: {overall_conclusion}")
    
 # 9. 改进建议
    print("\n9️⃣ 改进建议")
    print("-" * 30)
    
    suggestions = []
    
    if not results['linearity']:
        suggestions.append("• 考虑添加多项式项或使用非线性回归")
    
    if not results['independence']:
        suggestions.append("• 检查数据收集过程，考虑时间序列模型")
    
    if not results['homoscedasticity']:
        suggestions.append("• 使用稳健标准误或对因变量进行变换")
    
    if not results['normality']:
        suggestions.append("• 考虑Box-Cox变换或使用非参数方法")
    
    if not results['multicollinearity']:
        suggestions.append("• 移除高度相关的变量或使用正则化方法")
    
    if not results['outliers']:
        suggestions.append("• 仔细检查异常值，考虑稳健回归方法")
    
    if suggestions:
        print("📝 具体改进建议:")
        for suggestion in suggestions:
            print(suggestion)
    else:
        print("当前示例规则未触发改进建议；这不排除遗漏的设定问题或预测风险。")
    
    return results

# 对住房模型进行完整诊断
final_results = regression_diagnostic_checklist(
    housing_model, df_housing, 
    ['area', 'floor', 'age'], 'price'
)

⚠️ 注意事项与最佳实践

warning 🚨 **回归诊断注意事项**

诊断顺序很重要：

先检查线性关系和异常值
再检查残差的分布特性
最后检查多重共线性

常见误区：

过度依赖单一诊断指标
忽视业务背景和数据收集过程
盲目追求完美的统计假设
不考虑样本量对检验结果的影响

最佳实践：

结合图形诊断和统计检验
优先考虑实际意义而非统计显著性
记录所有诊断步骤和决策过程
对关键模型进行交叉验证

📚 扩展学习资源

理论深化

📖 《Applied Linear Regression》 - Weisberg
📖 《回归分析导论》 - 蒙哥马利等
🎓 在线课程：MIT《应用统计学》

实践工具

🐍 Python库：Statsmodels、Scikit-learn、Seaborn
📊 专业软件：R、SAS、SPSS
🔧 诊断包：influence_plot、diagnostic_plots

进阶主题

线性回归 - 回归诊断的基础
非线性回归 - 处理非线性关系
稳健回归 - 处理异常值的方法

tip 💡 **学习建议**

建立系统的诊断流程和检查清单
多练习真实数据的诊断分析
掌握图形诊断的解读技巧
理解各种违反假设的后果和处理方法
培养批判性思维，质疑模型结果

#统计学 #回归诊断 #模型检验 #假设检验 #异常值检测 #多重共线性 #Python统计分析