Python 逻辑回归：statsmodels 与 scikit-learn 实现

因变量是「良性/恶性」「存活/死亡」这类二分类结果时，线性回归直接失效：它输出连续的预测值，会给出 −0.3 或 1.7 这种不可能的概率，残差方差也随预测值变化，违反同方差假设。逻辑回归（logistic regression）先对线性组合做一次 logit 变换（logit transformation），把取值范围从整个实数轴压缩到 (0, 1)，再按最大似然（maximum likelihood）估计参数。

Python 里有两个入口。statsmodels 给出系数、标准误、置信区间，用于推断；scikit-learn 给出预测器和调参工具，用于预测。两者的系数含义一致，但默认行为有关键差别，最后一节单独说。

拟合模型

sklearn.datasets.load_breast_cancer 是威斯康星乳腺癌诊断数据，569 个样本、30 个特征，任务是判断肿瘤为恶性还是良性。取两个连续变量做演示：

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from sklearn.datasets import load_breast_cancer

df = load_breast_cancer(as_frame=True).frame
df = df[["mean radius", "mean texture", "target"]].copy()
df.columns = ["radius", "texture", "target"]

# sklearn 里 target = 0 是恶性、1 是良性，这里改写成 malignant = 1
df["malignant"] = 1 - df["target"]
df = df.drop(columns="target")

print(df.head(3))
print(df["malignant"].value_counts().sort_index())

   radius  texture  malignant
0   17.99    10.38          1
1   20.57    17.77          1
2   19.69    21.25          1
malignant
0    357
1    212
Name: count, dtype: int64

569 个样本里 212 个恶性、357 个良性，恶性占 37.3%。这个比例后面算准确率时要用到。

公式接口的写法和 R 的 glm(y ~ x, family = binomial) 对应，smf.logit() 已经把 family 固定成二项：

fit = smf.logit("malignant ~ radius", data=df).fit()
print(fit.summary())

Optimization terminated successfully.
         Current function value: 0.289992
         Iterations 8
                           Logit Regression Results
==============================================================================
Dep. Variable:              malignant   No. Observations:                  569
Model:                          Logit   Df Residuals:                      567
Method:                           MLE   Df Model:                            1
Date:                Wed, 23 Sep 2026   Pseudo R-squ.:                  0.5608
Time:                        07:19:39   Log-Likelihood:                -165.01
converged:                       True   LL-Null:                       -375.72
Covariance Type:            nonrobust   LLR p-value:                 1.192e-93
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept    -15.2459      1.325    -11.509      0.000     -17.842     -12.649
radius         1.0336      0.093     11.100      0.000       0.851       1.216
==============================================================================

前三行是优化器的运行信息，Date 和 Time 两行是本次运行的时间戳，都不参与统计结论。不想看到优化信息，写 .fit(disp=0)。

输出和 sm.OLS() 的 summary 结构接近，但有四处必须分清：

统计量列是 z 而不是 t。逻辑回归的参数靠迭代求解（Iterations 8 是迭代次数），没有 lm 那种精确的 t 分布，用的是标准正态近似。这意味着精度依赖样本量，n 很小的时候这一列的 p 值偏乐观。

系数是 log odds（对数优势），不是概率的变化量。radius 的系数 1.0336 表示「平均半径每增加 1，log odds 增加 1.0336」。这个数字无法直接讲给人听。

没有 R²。最大似然估计不存在「解释了多少方差」这个量，替代它的是偏差（deviance）。零模型（只含截距）的 LL-Null 是 −375.72，加入 radius 后 Log-Likelihood 升到 −165.01，Pseudo R-squ. 就是由这两者算出的 McFadden 伪 R²：

print(round(fit.prsquared, 4))
print(round(fit.llf, 3), round(fit.llnull, 3))

0.5608
-165.005 -375.72

0.5608 只能和同一批数据上的同类模型比较，不能说成「解释了 56.1% 的变异」，写论文时不要按 R² 的口径描述。AIC 用于模型比较，越小越好，第一个模型的 AIC 是 334.011。

LLR p-value 是似然比检验的 p 值，检验的是「当前模型是否优于零模型」。这里 1.192e-93，远小于 0.05。

把系数转成优势比

优势比（odds ratio，OR）是解读逻辑回归的唯一实用入口，取指数即可：

print(np.exp(fit.params))

Intercept    2.392227e-07
radius       2.811136e+00
dtype: float64

radius 的 OR 是 2.811。含义是：平均半径每增加 1 个单位，恶性的优势（odds）变为原来的 2.81 倍。系数为正对应 OR 大于 1，两者方向一致。

「每增加 1 个单位」这个说法不能省。半径的取值范围是 6.98 到 28.11，跨度和解释方式直接相关：换成每增加 2 个单位，OR 是 7.903；每增加 5 个单位，OR 涨到 175.55。同一个模型能给出三个差别巨大的数字，只报一个 OR 而不写清楚单位，读者无法判断效应大小。

置信区间同样要取指数：

print(np.exp(fit.conf_int()).round(3))

               0      1
Intercept  0.000  0.000
radius     2.342  3.374

radius 的 OR 区间是 [2.342, 3.374]，整段都在 1 的右侧，和 p 值给出的结论一致：半径与恶性概率正相关，方向明确。

截距的 OR 是 2.39e-07，置信区间上下限都被舍入成 0.000。这个数对应「半径为 0」时的恶性优势，而半径为 0 的肿瘤在数据里不存在。截距在这样的模型里通常没有实际含义，不必解读。

预测概率

fit.predict() 返回的是概率，不是 log odds。这一点和 R 的 predict() 不同——R 默认返回 log odds，必须显式写 type = "response"：

p = fit.predict(df)
print(p.head(8).round(3))

0    0.966
1    0.998
2    0.994
3    0.031
4    0.997
5    0.085
6    0.974
7    0.254
dtype: float64

前 3 个样本的半径都在 17 以上，模型给出 96% 以上的恶性概率；第 4 个样本半径小，概率只有 3.1%。概率随半径单调上升，符合系数为正的预期。

预测新样本时传入一个只含自变量列的 DataFrame，列名必须与拟合时一致：

new = pd.DataFrame({"radius": [10.0, 14.0, 20.0]})
print(fit.predict(new).round(4))

0    0.0073
1    0.3153
2    0.9956
dtype: float64

半径 10 的样本只有 0.7% 的恶性概率，半径 20 的样本是 99.6%。半径 14 落在中间，模型给出 31.5%——这类样本正是诊断中最难处理的部分，模型本身无法给出确定答案。

混淆矩阵与准确率

概率要变成分类结论，需要选一个阈值。最简单的是 0.5：

pred = (p > 0.5).astype(int)
print(pd.crosstab(pred, df["malignant"], rownames=["predicted"], colnames=["actual"]))
print("accuracy:", round((pred == df["malignant"]).mean(), 4))

actual       0    1
predicted
0          333   45
1           24  167
accuracy: 0.8787

569 个样本判对 500 个，准确率 87.87%。拆开看两个方向比单个准确率有用得多：

• 特异度（specificity）：357 个良性样本判对 333 个，93.28%
• 灵敏度（sensitivity，也叫召回率）：212 个恶性样本判对 167 个，78.77%
• 精确率（precision）：模型判为恶性的 191 个样本里 167 个确实是恶性，87.43%

两个错误方向的代价完全不同。45 个恶性样本被漏判成良性，在临床筛查里就是 45 个漏诊；24 个良性样本被误判成恶性，对应的是 24 次不必要的复查。选阈值时想清楚哪一类错误更不可接受，而不是默认 0.5。

statsmodels 也提供了现成的列联表，但行列方向与 pd.crosstab 相反：fit.pred_table() 的行是实际类别、列是预测类别。写论文附混淆矩阵时，方向标错会让灵敏度和特异度整个对调，两种写法选一种用到底。

阈值怎么选

0.5 不是定理，只是一个默认值。它假设两类同等重要，且模型输出的概率是校准过的。改变阈值对结果的影响：

for t in (0.3, 0.5, 0.7):
    pr = (p > t).astype(int)
    print(f"threshold {t:.1f}  accuracy {(pr == df['malignant']).mean():.4f}  flagged {pr.sum()}")

threshold 0.3  accuracy 0.8541  flagged 235
threshold 0.5  accuracy 0.8787  flagged 191
threshold 0.7  accuracy 0.8647  flagged 151

准确率在 0.5 处最高（0.8787），但这不是选 0.5 的理由。看另外两个数字：阈值降到 0.3，被判为恶性的样本从 191 增到 235，漏诊减少而误诊增加；阈值升到 0.7，判为恶性的只剩 151，误诊减少而漏诊增加。准确率对这两个方向的变化并不敏感，因为它把两类错误同等计数。

更有用的看法是精确率与召回率的权衡：

from sklearn.metrics import precision_recall_curve

prec, rec, thr = precision_recall_curve(df["malignant"], p)
print("PR 曲线上的点数:", len(prec))

PR 曲线上的点数: 457

precision_recall_curve 返回的 prec 和 rec 长度比 thr 多 1，最后一个是「全部判为负类」的边界点，画图时要注意对齐。

想同时看灵敏度与特异度的权衡，用 ROC 曲线（receiver operating characteristic curve）：

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from sklearn.metrics import roc_curve, roc_auc_score

fpr, tpr, thr = roc_curve(df["malignant"], p)
auc = roc_auc_score(df["malignant"], p)

fig, ax = plt.subplots(figsize=(5, 5))
ax.plot(fpr, tpr, lw=1.5, label=f"AUC = {auc:.3f}")
ax.plot([0, 1], [0, 1], color="grey", lw=1, ls="--")
ax.set_xlabel("1 - specificity")
ax.set_ylabel("sensitivity")
ax.legend(loc="lower right")
fig.savefig("roc.png", dpi=100, bbox_inches="tight")

图里横轴是 1 − 特异度，纵轴是灵敏度，虚线是随机猜测的对角线。曲线越靠近左上角，模型越好。

print("AUC:", round(auc, 4))

AUC: 0.9375

AUC 是曲线下面积，取值 0.5（无区分能力）到 1.0（完美区分）。0.9375 说明这个单变量模型已经能把两类样本较好地分开。AUC 的优势是不依赖阈值——它衡量的是模型给正类打高分的排序能力，因此拿它做模型间比较比拿准确率稳健得多。

加变量与模型比较

加入纹理均值再拟合一次：

fit2 = smf.logit("malignant ~ radius + texture", data=df).fit(disp=0)
print("Pseudo R-squ: %.4f    AIC: %.3f" % (fit2.prsquared, fit2.aic))
print("AIC fit1: %.3f" % fit.aic)

Pseudo R-squ: 0.6126    AIC: 297.123
AIC fit1: 334.011

伪 R² 从 0.5608 升到 0.6126，AIC 从 334.011 降到 297.123，两个指标都指向「加变量更好」。用似然比检验（likelihood ratio test）给出一个正式的 p 值：

from scipy import stats

lr = 2 * (fit2.llf - fit.llf)
print("logLik: %.3f vs %.3f" % (fit.llf, fit2.llf))
print("LR statistic %.3f  p %.3e" % (lr, stats.chi2.sf(lr, 1)))

logLik: -165.005 vs -145.562
LR statistic 38.888  p 4.489e-10

似然比统计量是两倍的对数似然差，在「只有截距的模型成立」这个原假设下服从自由度为 1 的卡方分布——自由度等于两个模型相差的参数个数。p = 4.489e-10，texture 的贡献显著。

print("corr: %.3f" % df["radius"].corr(df["texture"]))

corr: 0.324

两个自变量的相关系数是 0.324，共线性不严重，系数可以按「控制另一个变量不变」来解释。相关系数超过 0.8 时，系数会变得不稳定，符号甚至可能翻转，那时要改用 VIF 之类的指标判断，具体做法见 statsmodels 线性回归。

statsmodels 还是 scikit-learn

同一个模型用 scikit-learn 拟合，系数和 statsmodels 不完全一样：

from sklearn.linear_model import LogisticRegression

X, y = df[["radius"]], df["malignant"]

clf = LogisticRegression(max_iter=1000).fit(X, y)
print("sklearn default : coef %.4f  intercept %.4f" % (clf.coef_[0][0], clf.intercept_[0]))

clf2 = LogisticRegression(C=1e6, max_iter=1000).fit(X, y)
print("sklearn C=1e6   : coef %.4f  intercept %.4f" % (clf2.coef_[0][0], clf2.intercept_[0]))

print("statsmodels MLE : coef %.4f  intercept %.4f" % (fit.params["radius"], fit.params["Intercept"]))

sklearn default : coef 1.0252  intercept -15.1272
sklearn C=1e6   : coef 1.0340  intercept -15.2521
statsmodels MLE : coef 1.0336  intercept -15.2459

默认的 LogisticRegression() 给出 1.0252，比最大似然估计的 1.0336 小。原因是它默认带 L2 正则化，C=1.0 是正则强度的倒数，系数会被往 0 的方向收缩。把 C 调到 1e6，正则化几乎失效，系数变成 1.0340，与 statsmodels 的差距只剩求解器精度。

这决定了两者的分工：要系数、标准误、置信区间、p 值，用 statsmodels；要预测、交叉验证、网格搜索，用 scikit-learn。用 sklearn 的系数写论文汇报效应大小是错的，它已经被正则化改过。反过来，用 statsmodels 做批量预测和调参也不合适，它没有提供那套工具。

分工的实操方法见 scikit-learn 入门，那里讲了 Pipeline 与交叉验证。

两个必须知道的坑

完全分离（complete separation）。如果某个自变量能把两类样本彻底分开，最大似然估计不存在，系数会朝无穷发散。构造一组能完全分开的数据就能看到：

sep = pd.DataFrame({"x": [1, 2, 3, 4, 5, 6, 7, 8], "y": [0, 0, 0, 0, 1, 1, 1, 1]})
m = smf.logit("y ~ x", data=sep).fit(disp=0)
print(m.params.round(2))
print(m.bse.round(1))

PerfectSeparationWarning: Perfect separation or prediction detected,
parameter may not be identified

Intercept   -178.57
x             39.70
dtype: float64
Intercept    129841.5
x             29016.4
dtype: float64

8 个样本、一个普通自变量，系数却涨到 39.70，标准误 29016.4——比系数本身大三个数量级。这种数字没有任何统计意义，PerfectSeparationWarning 是唯一线索，而它只是一条警告，程序会照常跑完。看到系数绝对值异常大、标准误离谱，先怀疑分离问题，处理办法是减少自变量、增加样本，或改用带惩罚的估计（sklearn 的正则化在这里反而是优点）。

样本量。逻辑回归的经验法则是「每个自变量至少对应 10 到 20 个较少发生的那一类事件」。这里的较少类是 212 个恶性样本，理论上能支持十几个自变量；但如果某个子组的恶性样本只有二三十个，再往里加自变量就是在拟合噪声。

类别不平衡时准确率几乎没有意义。这份数据有 62.74% 是良性，一个把所有样本都判为良性的模型：

print("all-benign accuracy: %.4f" % (df["malignant"] == 0).mean())

all-benign accuracy: 0.6274

什么都不学就能拿到 62.74% 的准确率，而它的灵敏度是 0。恶性肿瘤一个都找不出来。数据不平衡更极端时这个数字还会更高，所以报告分类结果时至少同时给出灵敏度、精确率和 AUC。

下一步

因变量是计数而不是二分类时，statsmodels 提供了 sm.Poisson 和 sm.NegativeBinomial，公式接口的写法与 smf.logit 一致。三组以上的均值比较见 Python 方差分析；线性回归的完整流程和残差诊断见 statsmodels 线性回归。R 语言的 R逻辑回归 用 glm() 做同一件事，可以对照 summary() 输出的结构差异。