别只看R²了!用Python的statsmodels库做回归诊断,F检验和t检验到底怎么看?

当你用Python的 statsmodels 跑完一个线性回归,看到满屏的统计量时,是不是经常有种"每个字都认识但连起来就懵"的感觉?R²看起来不错,但F统计量、t统计量、p值这些到底在说什么?今天我们就来拆解这份"天书",让你真正看懂模型在说什么。

1. 回归诊断:从结果报告开始

运行 model.summary() 后,你会看到类似这样的输出(以statsmodels输出为例):

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.823
Model:                            OLS   Adj. R-squared:                  0.821
Method:                 Least Squares   F-statistic:                     487.5
Date:                Wed, 01 Jan 2020   Prob (F-statistic):           2.11e-74
Time:                        00:00:00   Log-Likelihood:                -1234.5
No. Observations:                 500   AIC:                             2475.
Df Residuals:                     497   BIC:                             2488.
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          2.4531      0.234     10.492      0.000       1.993       2.913
x1             0.4235      0.021     20.167      0.000       0.382       0.465
x2            -1.2345      0.112    -11.022      0.000      -1.455      -1.014
==============================================================================
Omnibus:                        2.345   Durbin-Watson:                   2.012
Prob(Omnibus):                  0.310   Jarque-Bera (JB):                2.456
Skew:                          -0.123   Prob(JB):                        0.293
Kurtosis:                       2.987   Cond. No.                         12.3
==============================================================================

1.1 关键指标速查表

统计量 位置 正常范围 解读要点
R-squared 顶部 0-1,越接近1越好 模型解释力,但容易受变量数量影响
Adj. R-squared 顶部 通常略低于R² 考虑变量数量后的修正R²
F-statistic 顶部 越大越好 模型整体显著性检验
Prob(F-statistic) 顶部 <0.05显著 模型整体p值
coef 中部 - 回归系数估计值
std err 中部 越小越好 系数估计的标准误差
t 中部 绝对值>2较显著 系数t检验统计量
P> t 中部

2. F检验:模型整体是否有效

F检验回答的核心问题是: 这个模型比只用均值预测强多少?

2.1 F统计量怎么算

F统计量的计算公式:

F = (解释的方差 / 预测变量数量) / (未解释的方差 / 自由度)
  = (ESS/df_model) / (RSS/df_resid)

在Python中,你可以这样验证:

# 计算F统计量的组成要素
ESS = np.sum((model.fittedvalues - y.mean())**2)  # 解释平方和
RSS = np.sum(model.resid**2)  # 残差平方和
TSS = ESS + RSS  # 总平方和

df_model = model.df_model  # 预测变量个数
df_resid = model.df_resid  # 残差自由度

F_statistic = (ESS/df_model)/(RSS/df_resid)
print(f"手动计算的F值: {F_statistic:.1f}, 模型报告的F值: {model.fvalue:.1f}")

2.2 如何解读F检验

  • Prob(F-statistic) :这是模型整体的p值
    • <0.05:模型至少有一个变量是显著的
    • 0.05:模型可能没有预测能力

注意:即使F检验显著,也不代表所有变量都重要,这就是为什么还需要t检验

3. t检验:每个变量是否重要

t检验关注的是: 在已有其他变量的情况下,这个变量还有必要存在吗?

3.1 t统计量的计算

每个系数的t值计算公式:

t = 系数估计值 / 标准误

在代码中验证:

# 获取第一个预测变量的结果
coef = model.params[1]  # 系数估计
std_err = model.bse[1]  # 标准误
t_value = coef / std_err

print(f"手动计算的t值: {t_value:.3f}, 模型报告的t值: {model.tvalues[1]:.3f}")

3.2 p值的实际含义

p值告诉你: 如果这个变量其实没用(真系数为0),看到当前结果的概率有多大?

  • p<0.05:有足够证据说明这个变量有用
  • p>0.05:不能证明这个变量有用(但不一定真的没用)

常见误区

  • p>0.05 ≠ 变量与y无关,只是在这个模型设置下不显著
  • p值受样本量影响极大,大样本时小效应也会显著

4. 实战案例:变量筛选策略

当你发现有些变量不显著时,该怎么办?来看一个真实案例:

4.1 向后剔除法步骤

  1. 拟合包含所有候选变量的完整模型
  2. 找出p值最大的变量
  3. 如果p>0.05,移除该变量
  4. 用剩余变量重新拟合模型
  5. 重复2-4直到所有变量都显著
import statsmodels.api as sm

def backward_elimination(data, target, predictors):
    """ 实现向后剔除法的函数 """
    while len(predictors) > 0:
        X = data[predictors]
        X = sm.add_constant(X)
        model = sm.OLS(target, X).fit()
        
        p_values = model.pvalues.drop('const')  # 排除截距项
        max_p = p_values.max()
        
        if max_p > 0.05:
            excluded = p_values.idxmax()
            predictors.remove(excluded)
            print(f"移除 {excluded} (p={max_p:.3f})")
        else:
            break
            
    return model

# 使用示例
import pandas as pd
from sklearn.datasets import load_diabetes

data = load_diabetes()
df = pd.DataFrame(data.data, columns=data.feature_names)
target = pd.Series(data.target)

final_model = backward_elimination(df, target, list(df.columns))
print("\n最终模型摘要:")
print(final_model.summary())

4.2 变量筛选注意事项

  • 不要一次性剔除多个变量 :变量之间可能有相关性,单独不显著的变量组合可能重要
  • 考虑理论重要性 :有些变量可能统计不显著但有理论价值,应该保留
  • 检查多重共线性 :高相关性的变量会互相"抢"显著性,VIF>10是警告信号

计算VIF的代码:

from statsmodels.stats.outliers_influence import variance_inflation_factor

def calculate_vif(X):
    """ 计算VIF的函数 """
    X = sm.add_constant(X)
    vif = pd.DataFrame()
    vif["Variable"] = X.columns
    vif["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
    return vif

print(calculate_vif(df))

5. 超越基础:高级诊断技巧

5.1 模型假设检查

线性回归有四大核心假设:

  1. 线性关系
  2. 误差项正态分布
  3. 同方差性
  4. 误差项独立

用Python检查这些假设:

import matplotlib.pyplot as plt
import seaborn as sns

# 1. 残差图检查线性性和同方差性
plt.figure(figsize=(10,6))
sns.residplot(x=model.fittedvalues, y=model.resid, lowess=True)
plt.xlabel("预测值")
plt.ylabel("残差")
plt.title("残差vs拟合值图")
plt.axhline(y=0, color='r', linestyle='--')

# 2. QQ图检查正态性
import scipy.stats as stats
plt.figure(figsize=(10,6))
stats.probplot(model.resid, plot=plt)
plt.title("残差QQ图")

# 3. 杜宾-沃森检验自相关
from statsmodels.stats.stattools import durbin_watson
print(f"杜宾-沃森统计量: {durbin_watson(model.resid):.2f}")
# 接近2表示无自相关,<1或>3有问题

5.2 离群值和强影响点

有些数据点会过度影响模型结果,需要识别:

# 计算影响度量
from statsmodels.stats.outliers_influence import OLSInfluence
influence = OLSInfluence(model)

# 库克距离 > 0.5的点需要关注
plt.figure(figsize=(10,6))
plt.stem(influence.cooks_distance[0], markerfmt=",")
plt.title("库克距离")
plt.xlabel("观测序号")
plt.ylabel("库克距离")

# 高杠杆点(hat values > 2*(p+1)/n 是阈值)
n, p = model.nobs, len(model.params)
high_leverage = influence.hat_matrix_diag > 2*(p)/n
print(f"高杠杆点数量: {sum(high_leverage)}")

6. 模型比较与选择

当你有多个候选模型时,如何选择最佳模型?

6.1 信息准则对比

准则 公式 特点
AIC -2LL + 2k 惩罚参数较少,适合预测
BIC -2LL + k*ln(n) 惩罚更重,适合解释
# 比较两个模型
model1 = sm.OLS(y, X1).fit()
model2 = sm.OLS(y, X2).fit()

print(f"模型1 AIC: {model1.aic:.1f}, BIC: {model1.bic:.1f}")
print(f"模型2 AIC: {model2.aic:.1f}, BIC: {model2.bic:.1f}")

6.2 交叉验证

更可靠的模型选择方法:

from sklearn.model_selection import cross_val_score
from sklearn.metrics import make_scorer
from sklearn.linear_model import LinearRegression

def rmse(y_true, y_pred):
    return np.sqrt(np.mean((y_true - y_pred)**2))

cv_scores = cross_val_score(LinearRegression(), X, y, 
                          cv=5, scoring=make_scorer(rmse))
print(f"RMSE交叉验证结果: {np.mean(cv_scores):.3f} ± {np.std(cv_scores):.3f}")

7. 常见问题解决方案

7.1 遇到多重共线性怎么办?

  • 移除相关性高的变量之一
  • 使用主成分分析(PCA)降维
  • 改用正则化方法(岭回归、Lasso)
# 岭回归示例
from sklearn.linear_model import Ridge
ridge = Ridge(alpha=1.0).fit(X, y)
print("岭回归系数:", ridge.coef_)

7.2 异方差性处理

  • 对因变量做变换(如对数变换)
  • 使用稳健标准误
  • 改用加权最小二乘法
# 使用稳健标准误
model_robust = sm.OLS(y, X).fit(cov_type='HC3')
print(model_robust.summary())  # 注意标准误的变化

7.3 非线性关系处理

  • 添加多项式项
  • 使用样条变换
  • 考虑广义加性模型(GAM)
# 多项式回归
from sklearn.preprocessing import PolynomialFeatures
poly = PolynomialFeatures(degree=2, include_bias=False)
X_poly = poly.fit_transform(X)
model_poly = sm.OLS(y, sm.add_constant(X_poly)).fit()

8. 完整分析流程示例

让我们用一个糖尿病数据集完整走一遍流程:

# 加载数据
from sklearn.datasets import load_diabetes
data = load_diabetes()
df = pd.DataFrame(data.data, columns=data.feature_names)
y = pd.Series(data.target)

# 1. 初始模型
X = sm.add_constant(df)
model = sm.OLS(y, X).fit()
print(model.summary())

# 2. 检查假设
# 残差图
plt.figure(figsize=(10,6))
sns.residplot(x=model.fittedvalues, y=model.resid, lowess=True)
plt.axhline(0, color='r', linestyle='--')

# 3. 变量筛选
final_model = backward_elimination(df.copy(), y, list(df.columns))

# 4. 检查最终模型
print("\n=== 最终模型 ===")
print(final_model.summary())

# 5. 模型诊断
influence = OLSInfluence(final_model)
plt.figure(figsize=(10,6))
plt.stem(influence.cooks_distance[0], markerfmt=",")
plt.title("库克距离 - 最终模型")

# 6. 预测新数据
X_new = final_model.model.exog[:5]  # 示例用前5个观测
predictions = final_model.get_prediction(X_new)
print("\n新数据预测结果:")
print(predictions.summary_frame())

记住,回归诊断不是一次性的工作,而是一个迭代过程。当发现问题时,修改模型后需要重新诊断,直到获得一个既简洁又有良好统计性质的模型。

更多推荐