statsmodels 线性回归：OLS 拟合与结果解读

Python 里做回归有两个常用库，分工不同。scikit-learn 面向预测，给你一组系数和预测值，不回答「这个系数显著吗」；statsmodels 面向统计推断，summary() 一张表里给出系数、标准误、t 值、p 值和置信区间。论文里要报的东西，statsmodels 都有。

本篇用 statsmodels.formula.api 的公式接口，写法与 R 的 lm(y ~ x) 基本一致，两边的输出数字也能对上。

准备数据

sklearn 自带的糖尿病数据集（diabetes）是标准的回归示例：442 位患者，10 个基线变量，目标变量是一年后疾病进展的量化指标。

import pandas as pd
from sklearn.datasets import load_diabetes
from sklearn.preprocessing import StandardScaler

raw = load_diabetes(as_frame=True).frame.copy()
raw.columns = [c.lower() for c in raw.columns]
print(raw.shape)
print(raw["target"].describe().round(2))

(442, 11)
count    442.00
mean     152.13
std       77.09
min       25.00
25%       87.00
50%      140.50
75%      211.50
max      346.00
Name: target, dtype: float64

这个数据集的自变量已经被 sklearn 做过 L2 归一化——每一列的平方和为 1，标准差只有 0.048。直接用原始列拟合，bmi 的系数会是 949 这种数字，因为它对应的是「归一化单位」而不是一个标准差，没法解读。

X = raw.drop(columns="target")
df = pd.DataFrame(
    StandardScaler().fit_transform(X), columns=X.columns
).assign(target=raw["target"].values)

print(df["bmi"].mean().round(6), df["bmi"].std().round(4))

-0.0 1.0011

StandardScaler 把每列变成均值 0、标准差 1 的 z 分数。这里打印出的 1.0011 不是误差：StandardScaler 按总体标准差（除以 n）归一化，而 pandas 的 .std() 默认按样本标准差（除以 n−1）计算，442 个样本下两者的比是 √(442/441) = 1.0011。

之后系数读作「自变量每增加 1 个标准差，因变量平均变化多少」，量纲统一，多个自变量的效应还能直接比大小。因变量保持原单位，这样预测值仍然是可以解释的分数。

从一个问题开始

先看体重指数（bmi）和疾病进展的关系：

import statsmodels.formula.api as smf

fit = smf.ols("target ~ bmi", data=df).fit()
print(fit.summary())

                            OLS Regression Results
==============================================================================
Dep. Variable:                 target   R-squared:                       0.344
Model:                            OLS   Adj. R-squared:                  0.342
Method:                 Least Squares   F-statistic:                     230.7
Date:                Wed, 23 Sep 2026   Prob (F-statistic):           3.47e-42
Time:                        07:23:41   Log-Likelihood:                -2454.0
No. Observations:                 442   AIC:                             4912.
Df Residuals:                     440   BIC:                             4920.
Df Model:                           1
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept    152.1335      2.974     51.162      0.000     146.289     157.978
bmi           45.1600      2.974     15.187      0.000      39.316      51.004
==============================================================================
Omnibus:                       11.674   Durbin-Watson:                   1.848
Prob(Omnibus):                  0.003   Jarque-Bera (JB):                7.310
Skew:                           0.156   Prob(JB):                       0.0259
Kurtosis:                       2.453   Cond. No.                         1.00
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

公式 target ~ bmi 读作「用 bmi 解释 target」。截距项默认包含；想去掉截距写 target ~ bmi - 1，除非有明确的物理理由，否则不要这么做，去掉截距会改变 R² 的定义，数值不再能和普通模型比较。

逐行读懂 summary 输出

这张表分四块，多数人只看中间的系数表。

左上角是模型的基本信息：因变量名、方法、样本量 No. Observations: 442、残差自由度 Df Residuals: 440（442 减 2 个待估参数）。Covariance Type: nonrobust 表示标准误没有做异方差稳健校正，后面单独说这一点。

coef 列是点估计。截距 152.13 表示 bmi 处于均值（z = 0）时，疾病进展指标的预测值就是 152.13——正好等于 target 的样本均值，因为自变量中心化后截距必然对应均值。bmi 的系数 45.16 是可解释的：bmi 每增加 1 个标准差，疾病进展指标平均上升 45.16 分。target 的标准差是 77.09，所以这个效应约为 0.59 个标准差。

std err 列是标准误，衡量这个系数估得有多准。t 列就是 coef / std err：45.16 / 2.974 = 15.19。P>|t| 列检验「该系数等于 0」，bmi 对应 0.000（实际是 3.5e-42 量级，表格只显示三位小数）。

最后一列 [0.025 0.975] 是 95% 置信区间，bmi 为 [39.32, 51.00]。区间的宽度比 p 值更有信息量：它给出效应的可能范围，而不是把结论简化成「有」或「没有」。单独取出来：

print(fit.conf_int().round(4))

                 0         1
Intercept  146.2894  157.9776
bmi         39.3159   51.0041

右上方是拟合优度：R-squared: 0.344、Adj. R-squared: 0.342、F-statistic: 230.7 配合 Prob (F-statistic): 3.47e-42。F 检验问的是「模型里所有系数是否全为 0」，单自变量时它和 bmi 的 t 检验等价（15.187² = 230.6）。AIC 和 BIC 用于比较模型，数值越小越好，但只在同一批数据上可比。

右下角是残差诊断量。Durbin-Watson: 1.848 接近 2，说明残差没有明显的自相关。Omnibus 和 Jarque-Bera 检验残差正态性，这里 p 值 0.003 和 0.026 都小于 0.05，残差正态性打了一点折扣——442 个样本下这个偏离不算严重，但如果要做精确的小样本推断就得留意。Cond. No. 是条件数，单变量模型恒为 1.00。

R² 与调整 R²

R² 的定义是「模型解释的方差占总方差的比例」。它有一个让人误用的特性：往模型里加任何变量，R² 都不会下降，哪怕新变量纯属噪声。

调整 R² 按自由度做了惩罚。加两个变量进去看变化：

fit2 = smf.ols("target ~ bmi + bp + s5", data=df).fit()
print(fit2.summary().tables[1])
print("R2", round(fit2.rsquared, 4), "Adj.R2", round(fit2.rsquared_adj, 4))

==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept    152.1335      2.653     57.342      0.000     146.919     157.348
bmi           28.6855      3.076      9.324      0.000      22.639      34.732
bp            12.4750      2.995      4.166      0.000       6.589      18.361
s5            25.8693      3.074      8.417      0.000      19.828      31.910
==============================================================================
R2 0.4801 Adj.R2 0.4765

summary().tables[1] 单独取出系数表，篇幅比整张 summary 短得多，贴进报告时更常这样用。R² 从 0.344 升到 0.480，调整 R² 从 0.342 升到 0.477，两者都在涨，说明新增的两个变量带来了实际解释力。

比较嵌套模型有一个更正式的做法，sm.stats.anova_lm() 给出 F 检验：

import statsmodels.api as sm

print(sm.stats.anova_lm(fit, fit2))

   df_resid           ssr  df_diff        ss_diff          F        Pr(>F)
0     440.0  1.719582e+06      0.0            NaN        NaN           NaN
1     438.0  1.362709e+06      2.0  356873.117068  57.352839  7.527665e-23

F = 57.35、p ≈ 7.5e-23，加入 bp 和 s5 确实显著改善了拟合。这一步要求两个模型嵌套（小模型的变量是大模型的子集）且用同一批数据拟合，否则 F 检验没有意义。fit2.compare_f_test(fit) 返回同样的 F 值和 p 值，可以直接用。

再看系数的变化：bmi 从单变量的 45.16 降到 28.69，几乎掉了一半。原因是 bmi 和 s5 的相关系数是 0.446，两者共享一部分解释力，单独放进模型时 bmi 把这部分也算在自己头上。这是多元回归里系数含义变化的典型来源——系数是「控制其他变量不变」的效应，不是单变量时的效应。

分类自变量：C() 与基准组

statsmodels 的公式接口不会自动把字符串列当分类变量，需要显式用 C() 包起来，否则报错。用鸢尾花数据集演示：

from sklearn.datasets import load_iris

iris = load_iris(as_frame=True).frame
iris.columns = ["sepal_len", "sepal_wid", "petal_len", "petal_wid", "species"]
iris["species"] = iris["species"].map({0: "setosa", 1: "versicolor", 2: "virginica"})

m_cat = smf.ols("sepal_len ~ C(species)", data=iris).fit()
print(m_cat.summary().tables[1])

============================================================================================
                               coef    std err          t      P>|t|      [0.025      0.975]
--------------------------------------------------------------------------------------------
Intercept                    5.0060      0.073     68.762      0.000       4.862       5.150
C(species)[T.versicolor]     0.9300      0.103      9.033      0.000       0.727       1.133
C(species)[T.virginica]      1.5820      0.103     15.366      0.000       1.379       1.785
============================================================================================

这是哑变量编码（dummy coding）：截距 5.006 是基准组 setosa 的均值，两个系数分别是 versicolor 和 virginica 相对 setosa 的均值差。三组均值就是 5.006、5.936、6.588。

基准组由水平的字母顺序决定，setosa 排在最前所以被选为基准。想换基准组，用 Treatment 显式指定：

print(smf.ols("sepal_len ~ C(species, Treatment(reference='virginica'))", data=iris)
      .fit().summary().tables[1])

==============================================================================================================================
                                                                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------------------------------------------------------
Intercept                                                      6.5880      0.073     90.492      0.000       6.444       6.732
C(species, Treatment(reference='virginica'))[T.setosa]        -1.5820      0.103    -15.366      0.000      -1.785      -1.379
C(species, Treatment(reference='virginica'))[T.versicolor]    -0.6520      0.103     -6.333      0.000      -0.855      -0.449
==============================================================================================================================

换基准组不改变模型拟合优度，三组均值仍然是 5.006、5.936、6.588，只是系数换了参照对象。注意行名带上了完整的 C(...) 表达式，打印出来很宽——把因子先转成有序的 pandas.Categorical，或者事后重命名结果表的行名，贴进报告时会好看很多。

有 k 个水平的因子会产生 k−1 个系数，剩下的那个「藏」在截距里。模型里加第二个分类变量时，statsmodels 默认同样是哑变量编码，不是方差分析式编码。

交互项：`x1:x2` 与 `x1*x2`

两个变量的效应不独立时需要交互项。* 同时展开主效应和交互项，: 只给交互项。

m_int = smf.ols("sepal_len ~ petal_len * C(species)", data=iris).fit()
print(m_int.summary().tables[1])
print("Adj.R2", round(m_int.rsquared_adj, 4))

======================================================================================================
                                         coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------------------------------
Intercept                              4.2132      0.407     10.341      0.000       3.408       5.018
C(species)[T.versicolor]              -1.8056      0.598     -3.017      0.003      -2.988      -0.623
C(species)[T.virginica]               -3.1535      0.634     -4.973      0.000      -4.407      -1.900
petal_len                              0.5423      0.277      1.959      0.052      -0.005       1.089
petal_len:C(species)[T.versicolor]     0.2860      0.295      0.969      0.334      -0.297       0.869
petal_len:C(species)[T.virginica]      0.4534      0.290      1.563      0.120      -0.120       1.027
======================================================================================================
Adj.R2 0.8349

每个系数都是相对 setosa 而言的。换算成三组各自的回归方程：

setosa：截距 4.2132，斜率 0.5423
versicolor：截距 4.2132 − 1.8056 = 2.4075，斜率 0.5423 + 0.2860 = 0.8283
virginica：截距 4.2132 − 3.1535 = 1.0597，斜率 0.5423 + 0.4534 = 0.9957

这正是分组各做一次回归的结果，交互模型只是额外假定了三组的残差方差相同。主效应不能单独解读：petal_len 的 0.5423 只是 setosa 一组的斜率，p = 0.052 说的是「setosa 组的斜率是否为零」，不能说成「花瓣长度整体不显著」。

报告交互模型时务必给出各组的分组斜率，或者画一张分组散点加拟合线。用 seaborn 的 lmplot 可以一行画出来，见 Seaborn统计可视化。

预测：置信区间与预测区间

new = pd.DataFrame({"bmi": [-1.0, 0.0, 1.0]})
print(fit.get_prediction(new).summary_frame(alpha=0.05).round(2).to_string())

     mean  mean_se  mean_ci_lower  mean_ci_upper  obs_ci_lower  obs_ci_upper
0  106.97     4.21          98.71         115.24        -16.17        230.12
1  152.13     2.97         146.29         157.98         29.13        275.14
2  197.29     4.21         189.03         205.56         74.15        320.44

mean 一列是点预测。后四列分成两组区间，宽度差很多：

mean_ci_lower / mean_ci_upper 是均值置信区间，描述回归线本身的不确定性。bmi 处于均值时是 [146.29, 157.98]，宽度约 11.7。
obs_ci_lower / obs_ci_upper 是观测预测区间，描述「某一位新患者」的不确定性。同一行是 [29.13, 275.14]，宽度约 246。

差别在于前者只包含参数估计误差，后者还要加上个体的随机误差。论文里报「模型预测的平均水平」用前者，报「某位患者的指标可能落在什么范围」必须用后者。把两者混用是审稿意见里常见的错误。

残差诊断

假设检验的前提是模型假设成立。下面用三变量的 fit2 检查残差，先看五数概括：

print(fit2.resid.describe().round(4).to_string())

count    442.0000
mean      -0.0000
std       55.5881
min     -140.2297
25%      -40.6387
50%       -2.1862
75%       38.2689
max      139.8018

中位数 −2.19 接近 0，残差大体居中；下四分位数 −40.64 与上四分位数 38.27 量级接近，没有明显偏斜。如果中位数明显偏离 0，说明有系统性偏差没被模型捕捉。

statsmodels 提供成套的诊断图：

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import statsmodels.graphics.regressionplots as rp

fig = rp.plot_regress_exog(fit2, "bmi", fig=plt.figure(figsize=(10, 8)))

plot_regress_exog 一次画出四张子图（左上为拟合线、右上为残差对拟合值、左下为偏回归图、右下为成分加残差图）。单独控制时可以分别用 sm.graphics.plot_fit(fit2, "bmi") 画拟合图，rp.plot_leverage_resid2(fit2) 画杠杆值对标准化残差。

看这几张图要找的东西与 R 一致：

残差对拟合值：理想情况是围绕 y = 0 的随机散点。出现 U 形说明线性假设不成立，需要加二次项；出现喇叭口（残差随拟合值变大而扩散）说明方差不齐，可以对因变量取对数，或者改用异方差稳健标准误。
Q-Q 图：点大致落在对角线上即可，两端轻微偏离可接受。严重弯曲说明残差非正态。
杠杆值对残差：看库克距离（Cook’s distance）大的点。杠杆高、残差又大的观测会明显拉偏回归线。

异常点别只凭肉眼判断，用数字更稳：

cooks = pd.Series(fit2.get_influence().cooks_distance[0], index=df.index)
print(cooks.sort_values(ascending=False).head(3).round(4).to_string())

256    0.0332
387    0.0239
289    0.0220

最大的库克距离只有 0.033，一般把大于 1 的点视为需要检查的对象，这份数据里没有。发现异常点不等于可以删掉它，先搞清楚它是不是录入错误；删点必须在方法部分写明。

前面 summary() 里显示的 Covariance Type: nonrobust 提醒标准误没做稳健校正。残差图出现喇叭口时，改用

fit2_robust = smf.ols("target ~ bmi + bp + s5", data=df).fit(cov_type="HC3")
print(fit2_robust.summary().tables[1])

==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept    152.1335      2.665     57.082      0.000     146.910     157.357
bmi           28.6855      3.279      8.749      0.000      22.259      35.112
bp            12.4750      3.197      3.902      0.000       6.208      18.742
s5            25.8693      3.081      8.396      0.000      19.830      31.908
==============================================================================

系数完全没变（异方差只影响标准误，不影响点估计），标准误全部变大：bmi 从 3.076 变成 3.279。表头也跟着变了——稳健协方差用的是渐近正态近似，所以统计量列从 t 改成 z、P>|t| 改成 P>|z|。正文里引用时要写清楚这一点，别一边用稳健标准误一边说是 t 检验。

多重共线性

自变量之间高度相关时，系数估计会变得不稳——标准误膨胀，甚至符号反转。数据集的 s1 和 s2 是两个相关系数 0.897 的血清指标，正好是现成的例子。先单独放 s1：

print(smf.ols("target ~ s1", data=df).fit().summary().tables[1])

==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept    152.1335      3.588     42.405      0.000     145.082     159.185
s1            16.3269      3.588      4.551      0.000       9.276      23.378
==============================================================================

再把 s2 一起放进去：

print(smf.ols("target ~ s1 + s2", data=df).fit().summary().tables[1])

==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept    152.1335      3.589     42.386      0.000     145.079     159.188
s1            21.9845      8.107      2.712      0.007       6.050      37.919
s2            -6.3096      8.107     -0.778      0.437     -22.244       9.625
==============================================================================

s1 的系数从 16.33 变成 21.98，标准误从 3.588 膨胀到 8.107，翻了一倍多；s2 的系数是负的且 p = 0.437，看起来「没有作用」。两个本来都显著的变量放进同一个模型后互相削弱，这是共线性的典型症状。

量化指标是方差膨胀因子（VIF，variance inflation factor）：

import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor

Xv = sm.add_constant(df[["s1", "s2"]])
for i, col in enumerate(Xv.columns):
    if col == "const":
        continue
    print(col, round(variance_inflation_factor(Xv.values, i), 3))

s1 5.102
s2 5.102

VIF 的含义是「该系数的方差被共线性放大了多少倍」。两个变量时 VIF 相同，都等于 1/(1 − 0.897²) = 5.10，与相关系数一一对应。经验上 VIF 大于 5 或 10 就该处理：

VIF	判断
< 5	基本不用担心
5 – 10	系数不稳，报结果时要说明
> 10	严重共线，必须处理

处理办法有三种：只留物理意义更明确的那个变量、把两者合并成一个综合指标、或者用主成分回归。不要因为共线性就随手删变量——如果两个变量都是研究问题需要的，删掉会引入遗漏变量偏误。共线性只影响系数的精度，不影响模型的预测能力。

与 scikit-learn 的分工

两个库对同一份数据给出完全相同的系数：

from sklearn.linear_model import LinearRegression

cols = ["bmi", "bp", "s5"]
sk = LinearRegression().fit(df[cols], df["target"])
print("statsmodels:", fit2.params[1:].round(4).tolist())
print("sklearn    :", sk.coef_.round(4).tolist())

statsmodels: [28.6855, 12.475, 25.8693]
sklearn    : [28.6855, 12.475, 25.8693]

两者都解最小二乘问题，所以点估计一致。差别在别处：

需求	用哪个
p 值、置信区间、F 检验、残差诊断	statsmodels
预测、交叉验证、正则化（Ridge / Lasso）	scikit-learn
与预处理步骤串联成流水线	scikit-learn

scikit-learn 的 cross_val_score 给出模型在未见数据上的表现，这是 statsmodels 不直接提供的能力：

from sklearn.model_selection import KFold, cross_val_score

cv = KFold(n_splits=5, shuffle=True, random_state=0)
scores = cross_val_score(sk, df[cols], df["target"], cv=cv, scoring="r2")
print(scores.round(4).tolist())
print("mean %.4f  sd %.4f" % (scores.mean(), scores.std()))

[0.3335, 0.4216, 0.4991, 0.4795, 0.5847]
mean 0.4637  sd 0.0835

五折的 R² 平均值 0.4637，比训练集上的 0.4801 低一些，这个差额就是过拟合的量。折与折之间在 0.334 到 0.585 之间波动，说明 442 个样本下模型表现对数据划分比较敏感。论文里报交叉验证结果时必须给出标准差，只报平均值会掩盖这种波动。

下一步

线性回归的因变量必须是连续变量。因变量是「是/否」这类二分类结果时，要换成逻辑回归，见 Python逻辑回归；三组以上的均值比较见 Python方差分析。R 语言的同一套回归逻辑在 R线性回归里，系数与 p 值和本篇完全一致，可以对照着看两种写法。