别只看R²了!用Python的statsmodels库做回归诊断,F检验和t检验到底怎么看?
·
别只看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 向后剔除法步骤
- 拟合包含所有候选变量的完整模型
- 找出p值最大的变量
- 如果p>0.05,移除该变量
- 用剩余变量重新拟合模型
- 重复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 模型假设检查
线性回归有四大核心假设:
- 线性关系
- 误差项正态分布
- 同方差性
- 误差项独立
用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())
记住,回归诊断不是一次性的工作,而是一个迭代过程。当发现问题时,修改模型后需要重新诊断,直到获得一个既简洁又有良好统计性质的模型。
更多推荐


所有评论(0)