1. 先想清楚：什么问题该交给logistic回归
1.1 从"是/否"的结局说起
我在实际项目里被问得最多的一类问题，结局都是一句话就能回答的：这个用户下个月会不会流失、这笔订单会不会被退回、这台设备在保期内会不会报故障、这个病人复查指标会不会转阳。结局只有两类，且我们真正关心的不是一个确定的0或1，而是"发生这件事的概率有多大"。这就是logistic回归最典型的战场。
R语言在这件事上的位置很特殊。它不像某些可视化软件那样点几下就出结果，但只要你愿意写几行代码，从数据清洗、模型拟合、系数解释到后续验证，整条链路都在同一个环境里完成，中间不需要导来导去。我这些年做二分类预测，凡是需要把"为什么这么判"讲清楚的场合，第一选择基本都是R里的glm。
这里要先划清一个边界：logistic回归不是"分类算法"，它是一个概率模型。模型输出的是事件发生的概率，分类只是你在概率上加了一道阈值切出来的副产品。搞清楚这一点，后面很多争论（比如"阈值取0.5还是0.3"）才有讨论的基础。很多人一上来就问准确率多少，其实方向就偏了。
这篇文章适合两类人：一是刚学完线性回归、准备接触二分类建模的朋友，跟着走一遍能建立正确的直觉；二是已经在用R跑模型、但每次看summary()输出都只盯着星号看的人——后面几节专门讲怎么把那些数字翻译成业务能听懂的话。数据、代码、踩过的坑，我都尽量摊开写。
提示：这篇文章定位是系列的第一篇，重点放在"数据怎么准备、模型怎么跑起来、系数怎么解释"这三件事上。模型验证的更细部分会在后面单独展开，这里只做最基础的检查，够你把第一版模型交出去。
1.2 为什么不能直接用线性回归硬套
很多人第一反应是：结局是0/1，那我把它当连续变量，跑一个lm()不就完了？技术上确实跑得出来，但出来的东西是错的，而且错得很隐蔽。
第一，线性回归的预测值没有边界。你拟合出来的直线可以给出-0.3或者1.8这样的预测值，而概率的取值范围只能是0到1。落在边界外的预测值没有任何解释意义，你没法跟人说"这个用户流失的概率是负的百分之三十"。
第二，线性回归默认残差方差恒定（同方差性）。但结局是0/1的时候，方差完全由均值决定：p越接近0或1，方差p(1-p)越小；p接近0.5时方差最大。这种方差随均值变化的结构直接违反了同方差假设，标准误的估计会失真，显著性检验也就不可信了。
第三，最关键的一点，logistic回归做的事情不是"拟合一条线"，而是对几率做对数变换，把(0,1)这个区间"拉直"成整条实数轴。直观理解：概率从0.01变到0.02，看起来只涨了1个百分点，但几率从0.0101涨到0.0204，翻了一倍；从0.50变到0.51，几率只从1.0变到1.04。logit变换正是把这种"两端变化剧烈、中间平缓"的特性抹平，让线性关系能在变换后的尺度上成立。
我把三种思路的差别整理成一张表，选型的时候对照看比较省事：
| 做法 | 预测值范围 | 系数解释 | 适用场合 |
|---|---|---|---|
lm(y ~ x)硬套 | 无界，可能超出[0,1] | 每单位x带来的概率绝对变化（常数） | 不推荐用于二分类结局 |
glm(..., family = binomial) | 概率(0,1)，logit尺度上无界 | 每单位x带来的对数几率变化，可转OR | 二分类结局的标准做法 |
| 判别分析 / 树模型 | 视方法而定 | 不同 | 分布假设不满足或关系高度非线性时考虑 |
还有一类折中方案是"线性概率模型 + 稳健标准误"，在计量经济学里确实有人这么干，系数可以直接读作边际效应，好解释。但它的预测值照样会越界，而且函数形式假设太强。我在实际项目里只在极少数需要快速看方向性结论的场合用过，最终交付还是回到logistic回归。
1.3 样本量与事件数：EPV这条隐形红线
这一条是新手最容易忽略、后果又最严重的。决定logistic回归能不能稳的，不是总样本量，而是较少那一类的事件数。
业内流传的经验法则是EPV（events per variable，每个自变量对应的事件数）至少10。也就是说，如果你的结局里"阳性"只有40例，那你最多放4个自变量进去。放多了会怎样？系数估计的偏差变大、标准误被低估、置信区间失真，运气不好的时候还会碰上完全分离（后面第6节专门讲）。模型在训练集上看着还行，换一批数据就崩。
举一个我实际遇到的场景：某团队做设备故障预测，总样本两万多条，听起来很充裕，但真正故障的只有63条。他们一口气放了二十多个特征，模型跑出来几个变量的OR值大到离谱（几十甚至上百），标准误也大得离谱。这不是模型发现了强关联，这是样本撑不住这么多参数。后来砍到5个变量，用Firth惩罚似然重跑，系数才回到合理区间。
样本量的估算可以用pwr包做一个粗略的规划：
library(pwr) # 双侧检验，显著性水平0.05，检验效能0.8 # h为效应量，二分类场景常用0.3~0.5之间 pwr.p.test(h = 0.3, sig.level = 0.05, power = 0.8, alternative = "two.sided")更贴近回归场景的可以用pwr里的pwr.f2.test近似，或者直接用R包WebPower、powerMediation做更细的规划。我的习惯是：先按EPV反推一个上限，再拿统计效能正着验算一遍，两个数取小的那个。宁可少放变量，也不要为了"模型看起来丰富"硬塞。
注意：EPV=10是经验值不是铁律。如果变量间相关性很低、事件数虽然少但分布均衡，EPV降到5~8有时也能接受；反过来，如果存在强共线性或者类别极度不平衡，EPV=10都未必够。这个数字的作用是提醒你别越界，不是给你一个精确的门槛。
2. 建模前的数据体检
2.1 结局变量的编码与类别不平衡
R里glm对结局的编码有讲究，这一点我见过太多人栽跟头。结局如果是因子（factor），R会把第一个水平当作"失败"（0），第二个水平当作"成功"（1）。而模型建模的是P(y = 第二个水平)。如果你的因子水平是按字母顺序排的，"no"在前"yes"在后，那正好；但如果水平是"正常"和"异常"，按拼音排序可能把"异常"排在前面，你算出来的OR就全部反过来了，而且不会报任何错。
稳妥的做法是显式指定：
df$y <- factor(df$y, levels = c("no", "yes")) # 第一个是无事件，第二个是事件 # 或者直接用0/1数值型，R会当成0=失败、1=成功 df$y <- as.integer(df$y == "yes")用0/1数值型反而最不容易出错，我在交付脚本里一般统一转成0/1，注释里写清楚1代表什么。
接下来是类别不平衡。真实业务里阳性率5%甚至1%以下很常见。不平衡本身不违反logistic回归的假设，极大似然估计在理论上仍然一致，但会带来两个实际麻烦：一是少数类的系数标准误很大，结论不稳；二是默认0.5阈值下的分类结果几乎全判成多数类，准确率能到95%以上却毫无用处。
我的处理顺序是这样的：先看不平衡到什么程度。阳性率在20%以上，基本不用管，正常建模即可；10%~20%，解释时留意一下置信区间；低于10%，就要考虑后续的阈值调整、重抽样或者加权。注意，不要在建模前无脑做SMOTE之类的过采样然后直接拟合，那会人为放大少数类，让预测概率整体偏高，概率校准就废了。如果确实要重抽样，放在交叉验证的每一折内部做，或者干脆用权重法。
# 用权重处理不平衡，比过采样更温和，且不改变概率尺度太多 w <- ifelse(df$y == 1, 1, sum(df$y == 1) / sum(df$y == 0)) fit_w <- glm(y ~ ., data = df, family = binomial, weights = w)权重的思路是给少数类更高的权重，让似然函数"更在意"它们。缺点是标准误的解释会变得微妙，严格来说不再是原始总体的标准误。我一般只在探索阶段用，最终交付模型还是用原始数据拟合，靠阈值和后处理来平衡。
2.2 缺失值、异常值与连续变量的形态
缺失值这件事，glm的默认行为是整行删除（na.omit）。这意味着你有10个变量，每个变量各缺5%，最后可能只剩六成样本。样本量一缩，EPV那条红线就又压过来了。
先看缺失的分布：
library(naniar) miss_var_summary(df) gg_miss_var(df)关键不是缺失比例，而是缺失是否与结局相关。如果"收入"这个字段缺失的人恰好是违约率特别高的一群人，那整行删除就是有偏的，删得越多偏得越狠。这种情况我一般会先加一个"是否缺失"的指示变量进去看看它显不显著，显著就说明缺失本身有信息量。真正要填补的话，mice包的多重插补是标准做法，流程是插补、分别建模、按Rubin规则合并结果：
library(mice) imp <- mice(df, m = 5, method = "pmm", seed = 2024) fit_imp <- with(imp, glm(y ~ age + sex + bmi, family = binomial)) pool(fit_imp)多重插补比单一插补（比如均值填充）好在它保留了填补的不确定性，标准误会更诚实。缺点是要跑五遍、流程复杂，团队里不熟悉的人接手容易出错。我的取舍是：缺失比例低于3%且随机，单次中位数插补加上缺失指示变量就够了；超过10%或者明显非随机，老老实实上mice。
连续变量的形态也要看。logistic回归假设的是logit尺度上的线性关系，不是原始尺度。这一点极容易被忽视：一个变量在原始尺度上看着跟结局没啥关系，取个分组或者做个样条之后关系就出来了。我常用的检查方式是把连续变量分箱（比如五分位），算每箱的事件率，再画出来看形状：
library(dplyr) df %>% mutate(bin = ntile(age, 5)) %>% group_by(bin) %>% summarise(n = n(), rate = mean(y == 1), mean_age = mean(age))如果事件率随分箱单调上升但两端变缓，经典的S形，那原始尺度进模型没问题；如果中间高两头低，或者出现明显的折点，就该考虑加二次项、分段或者用splines::ns()做限制性立方样条。我个人的经验是，年龄、收入、时长这类变量几乎从来不是纯线性的，随手加个样条往往能带来实打实的提升，而且解释起来也不算太难——用predict画一条预测概率曲线给业务方看，比解释系数直观得多。
异常值方面，连续变量里那种极端大值对logit尺度的影响比线性回归更剧烈，因为它在logit上可能把预测概率推到接近0或1。我的做法是先算标准化残差和Cook距离，把影响点标出来，再逐个判断是录入错误还是真实极端值：
plot(fit, which = 4) # Cook距离 plot(fit, which = 5) # 标准化残差 vs 杠杆值真实极端值不要随手删，删了就是人为制造乐观结果。可以做一个敏感性分析：删掉之后再跑一遍，看系数变化大不大，变化大就在报告里说明。
2.3 共线性与变量筛选的取舍
共线性在logistic回归里的表现和线性回归类似，但后果更容易被误读：两个高度相关的变量同时进模型，各自的系数可能都不显著，甚至符号反了，可它们合起来是有预测力的。新手看到"都不显著"就把两个都删了，等于把有用的信息一起扔掉。
先量化：
library(car) vif(fit)VIF超过5要留意，超过10基本就得处理了。但要注意，VIF对二分类变量和只有两三个水平的因子不太适用，别看到因子变量的VIF数值很奇怪就慌。
处理共线性有三条路，我按使用频率排：
第一，合并。两个变量测的其实是同一件事（比如"月收入"和"月消费"），那就选业务上更稳定、缺失更少的那个，或者做个比值/差值构造新变量。这比统计筛选靠谱得多，因为背后有业务逻辑支撑。
第二，降维。变量多且高度相关时，主成分或者偏最小二乘能压掉共线性，但代价是系数没法解释了。我在需要"讲清楚每个因素的作用"的项目里不用这招，只在纯预测场景下考虑。
第三，正则化。glmnet的弹性网络在共线性严重、变量又多的场合非常好用，L2项能把相关变量的系数往一起拉，L1项还能顺便做筛选：
library(glmnet) x <- model.matrix(y ~ . - 1, data = df) y <- df$y cv_fit <- cv.glmnet(x, y, family = "binomial", alpha = 0.5) plot(cv_fit) coef(cv_fit, s = "lambda.1se")这里alpha在0和1之间，0是岭回归，1是lasso，0.5左右是弹性网络，兼顾筛选和稳定。我用下来，lambda.1se比lambda.min更保守，留下来的变量少但更稳，适合交付。
关于变量筛选，我要多说一句反直觉的话：不要用单因素筛选（逐个跑单变量logistic，p<0.05才进多因素）作为标准流程。这个做法在流行病学里流传很广，但它有两个硬伤：一是会漏掉那些单独看没作用、但控制其他变量后才显现的抑制变量；二是多重检验没有校正，假阳性率高。更合理的做法是基于业务逻辑和领域知识先定候选集，再用AIC/BIC或者交叉验证做后向筛选，同时把"变量数不超过EPV红线"作为硬约束。
| 筛选方式 | 优点 | 风险 | 我的使用建议 |
|---|---|---|---|
| 业务先行 + 后向AIC | 可解释、稳定 | 依赖领域知识 | 首选，绝大多数项目用这个 |
| 单因素p值筛选 | 操作简单 | 漏变量、假阳性 | 只作为初筛参考，不作为最终依据 |
| Lasso/弹性网络 | 自动处理高维共线 | 系数有偏，需去偏修正 | 变量多于事件数十分之一时用 |
| 全模型不筛选 | 无选择偏倚 | 过拟合、EPV超标 | 仅在变量数极少（3~5个）时 |
3. 在R里把模型跑起来
3.1 glm函数的关键参数
核心就一行，但每个参数都值得说清楚：
fit <- glm(y ~ age + sex + bmi + smoke, data = df, family = binomial(link = "logit")) summary(fit)family = binomial指定了二项分布，link = "logit"是连接函数。其实写family = binomial默认就是logit链接，但我在交付脚本里习惯写全，别人接手时一眼能看明白。另外两个链接函数偶尔会用到：probit在计量经济学里多，cloglog适合处理"一段时间内是否发生"的场景，解释上略有不同，这里不展开。
data参数指向数据框，公式里的变量名会自动在这个数据框里找。不要在公式里写df$age这种形式，会造成后续predict时的环境混乱——这是我早期踩过的坑，模型对象里保存的公式带上了数据框名字，换一份数据做预测时报错报得莫名其妙。
其他几个参数在实际项目里会用到：
weights用于前面说的不平衡处理，或者当你的数据是汇总形式时（每个观测代表n个人，其中k个事件），用来传入试验次数。
offset用于处理不同观测的"暴露量"不同。比如预测某段道路是否发生事故，不同路段的车流量差几十倍，把log(车流量)作为offset放进去，模型就变成"每单位暴露的事件率"，而不是绝对事件数。这个参数用的人少，但在比率型结局里非常关键。
subset和na.action我一般不用，数据准备阶段就处理干净，模型里保持简洁。倒是control参数偶尔要调，默认的最大迭代次数25次在数据尺度差异很大时可能不够：
fit <- glm(y ~ ., data = df, family = binomial, control = glm.control(maxit = 50, epsilon = 1e-8))变量尺度差异大（比如一个变量取值0~1，另一个0~100000）会引起收敛慢甚至数值问题，标准的处理是先把连续变量标准化。scale()一行搞定，而且标准化之后的系数大小可以直接横向比较重要性——虽然严格来说还要考虑变量的实际取值范围，但比原始尺度直观多了。
df_z <- df df_z[c("age", "bmi", "income")] <- scale(df[c("age", "bmi", "income")])注意标准化之后再解释OR，单位就变成"每增加一个标准差"了，写报告的时候一定要标注清楚，否则读者会按"每增加1岁"去理解，差之毫厘谬以千里。
3.2 从summary里读什么
summary(fit)的输出信息量很大，但新手常犯的错是只看最后那列p值。我按重要性给你排个序。
第一看收敛状态和残差偏差。输出最后如果有"Number of Fisher Scoring iterations"，说明正常收敛了；如果出现"algorithm did not converge"或者"fitted probabilities numerically 0 or 1 occurred"，别急着往下读，先解决收敛问题（第6节）。
第二看系数表。四列分别是估计值Estimate、标准误Std. Error、z值、p值。这里要建立一个概念：Estimate是在对数几率（log-odds）尺度上的，不是概率尺度，也不是OR。z值就是Estimate除以标准误，检验的原假设是"这个系数等于0"，在logit尺度上等于说"这个变量对应的OR等于1"。
第三看零偏差和残差偏差。Null deviance是只用截距项时的偏差，Residual deviance是加了变量之后的偏差。两者之差近似服从卡方分布，自由度是自变量个数，可以做一个整体的似然比检验：
# 整体模型检验 1 - pchisq(fit$null.deviance - fit$deviance, df = fit$df.null - fit$df.residual)这个检验回答的是"这组变量作为一个整体有没有用"，比逐个看p值更靠谱。至于模型之间比较，用anova(fit1, fit2, test = "Chisq")做偏差分析，比看AIC更严谨（AIC适合非嵌套模型比较，嵌套模型用似然比检验更合适）。
第四看AIC。同一份数据、同一个结局下，AIC越小越好，但只能横向比，绝对值没有意义。我一般用它做后向筛选的准则，step(fit, direction = "both")，但注意step是按AIC走的，会倾向于保留较多变量，跑完之后要人工再看一遍，把业务上说不通的变量拿掉。
第五，summary不显示的东西才是重点：AUC、校准、各变量的实际效应大小。这些后面几节讲。
提示：p值大不等于变量没用。在样本量小的项目里，一个OR=2.5但p=0.08的变量，很可能只是效能不足，而不是真的没关联。看置信区间比看p值信息量大得多，如果区间宽得离谱，说明是数据不够，不是效应不存在。
3.3 预测概率与阈值设定
模型跑通之后，第一步是把线性预测值转成概率。predict的type参数有三个常用取值，搞混了就会得到莫名其妙的结果：
# 线性预测值（log-odds尺度） lp <- predict(fit, newdata = df, type = "link") # 预测概率，最常用 phat <- predict(fit, newdata = df, type = "response") # 直接预测类别，按0.5切 cls <- predict(fit, newdata = df, type = "class")type = "response"做的事情就是把lp过一遍exp(x)/(1+exp(x))。手动验证一下有助于建立直觉：如果某人的线性预测值是0，概率就是0.5；是1，概率约0.731；是-1，概率约0.269；是2，概率约0.881。线性预测值每增加1，几率变成原来的e≈2.718倍——这句话是整个logistic回归解释体系的基石，记住它后面就顺了。
type = "class"我基本不用，因为它固定用0.5切，而0.5几乎从来不是业务上的最优点。阈值该取多少，取决于两类错误的代价。举个虚构的例子：某团队做欺诈识别，漏判一笔的损失远大于误判一笔带来的客服成本，那阈值就该往下压，比如0.2，宁可多拦一些再人工复核。
选阈值的常用方法是看ROC曲线上离左上角最近的点（约登指数最大），或者直接按业务能承受的人工复核量倒推：
library(pROC) r <- roc(df$y, phat, quiet = TRUE) best <- coords(r, "best", ret = c("threshold", "sensitivity", "specificity")) best更贴近业务的做法是算一下不同阈值下的预期成本，取最小的那个。这个思路我在实际项目里用得最多，因为它把"模型指标"翻译成了"钱"，决策者一听就懂。
| 阈值 | 灵敏度 | 特异度 | 预计人工复核量 | 适用场景 |
|---|---|---|---|---|
| 0.5 | 低 | 高 | 少 | 后果对称、追求省事 |
| 0.3 | 中 | 中 | 中 | 漏判代价略高于误判 |
| 0.15 | 高 | 低 | 多 | 漏判代价极高、有人工兜底 |
还有一个容易被忽略的点：用训练集数据选出来的阈值是乐观的。正确做法是在交叉验证或者独立验证集上选阈值，再拿到测试集上评估。如果数据量实在不够，至少要做一次bootstrap重抽样看阈值的波动范围。
4. 系数怎么翻译成人话
4.1 log-odds、OR与置信区间
这一步是logistic回归真正的门槛。系数本身没法直接读，必须做指数变换变成OR（比值比）。完整的一套代码：
# 系数、OR、95%置信区间 res <- cbind( beta = coef(fit), OR = exp(coef(fit)), CI_low = exp(confint.default(fit)[, 1]), CI_high = exp(confint.default(fit)[, 2]) ) round(res, 3)confint.default用的是Wald法（估计值±1.96×标准误再取指数），快但对小样本或者系数较大时不准；confint用的是剖面似然法，更准但慢一些。样本量小的时候我倾向用后者：
exp(confint(fit)) # 剖面似然置信区间，更可靠OR怎么读？一句话：在其他变量保持不变的前提下，该变量每增加一个单位，事件发生的几率变成原来的OR倍。OR=1.5意味着几率增加50%，OR=0.6意味着几率降低40%。注意是"几率"不是"概率"，这个区别在结局发生率较高（比如超过10%）时非常明显。OR=2在基线概率10%时，对应的概率从10%涨到约18%；在基线概率1%时，只涨到约2%。所以报告OR的时候，最好同时给出基线风险，否则读者会高估效应。
置信区间比p值重要。区间不跨1就等价于p<0.05（Wald法下严格等价），但区间还能告诉你效应的可能范围。一个OR=1.8、区间[1.05, 3.1]的变量，和一个OR=1.8、区间[1.4, 2.3]的变量，虽然点估计一样，可信程度完全不同。我在给业务方汇报时，通常直接画一张森林图，把OR和区间一起摆出来：
library(sjPlot) plot_model(fit, show.values = TRUE, value.offset = 0.3, transform = "exp")这张图比一张系数表好用十倍，尤其在有十几个变量的时候，谁强谁弱一眼就看出来了。
4.2 连续变量的OR不是"每增加1个单位"那么简单
这是我最常被追问的点。"年龄的OR是1.03，那是不是说年龄每大一岁，风险涨3%？"——数学上没错，但这么讲几乎等于没讲，因为1.03这个数字太接近1，听众会直觉觉得"没啥影响"，而实际上从30岁到60岁，累积效应是1.03的30次方，约2.43倍。
所以对连续变量，我一般会换算成更有业务意义的单位：年龄按10岁一档、收入按万元一档、时长按小时一档。做法很简单，把变量先除以单位再进模型，或者事后直接算OR^k：
# 年龄每增加10岁的OR OR_10 <- exp(coef(fit)["age"] * 10) # 或者用模型直接算 exp(coef(fit)["age"] * 10)除了OR，还有一种更直观的表达叫平均边际效应（AME），意思是"在其他变量取实际值的情况下，这个变量增加一单位，概率平均变化多少个百分点"。margins包一行搞定：
library(margins) m <- margins(fit) summary(m)AME的单位是概率，业务方一听就懂："年龄每增加一岁，流失概率平均上升0.4个百分点"。缺点是它依赖于样本中其他变量的分布，换一个人群结论会变；而且对同一个变量，AME在不同基线风险下并不恒定。我的习惯是两种都给：OR用于说明方向和相对强度，AME用于说明绝对量级，配合着讲最清楚。
4.3 交互项与分层解读
交互项在logistic回归里尤其容易翻车，因为交互项的系数不是OR，是指数化之后得到的"OR的比值"（ratio of odds ratios）。这一点跟线性回归里交互项可以相对直观地理解为"斜率之差"完全不同。
举例，模型里有sex和smoke的交互项。指数化之后得到1.8，正确读法是：吸烟对结局的OR，在两组性别之间的比值是1.8，也就是某一组的OR是另一组的1.8倍。不要读成"吸烟且某性别的OR是1.8"。这个错误我在审阅别人的报告时见过不止一次。
正确展示交互效应的方法是按其中一个变量分层算简单效应：
library(emmeans) emm <- emmeans(fit, ~ smoke | sex, type = "response") emm # 或者在对数几率尺度上看差异 pairs(emmeans(fit, ~ smoke | sex))还有一个技术细节：如果模型里有交互项，那么主效应的系数含义变了。smoke的系数不再代表"吸烟的整体效应"，而是"在sex取参照水平时的效应"。所以要么把交互项去掉，要么在报告里明确说明参照组是什么，否则读者会误解。
连续变量之间的交互更难处理，因为乘积项和主效应高度相关，会导致VIF飙高。标准做法是先中心化：
df$age_c <- df$age - mean(df$age) df$bmi_c <- df$bmi - mean(df$bmi) fit_int <- glm(y ~ age_c * bmi_c + sex, data = df, family = binomial)中心化不改变交互项的系数，但让主效应变得可解释（在另一个变量取均值时的效应），还能缓解一部分共线性。这个技巧在实操中非常有用，代价只是解释时要多说一句"其他变量取均值时"。
我的经验是：交互项要克制。每加一个交互项，解释复杂度上升一个台阶，而且很容易过拟合。除非有明确的理论假设或者数据里能看到清晰的分层模式，否则不要在探索阶段乱试一堆交互。真要试，就用anova()比较加与不加的偏差差异，显著了再留。
5. 模型好不好用：评估与验证
5.1 区分度：ROC与AUC
AUC衡量的是"随机抽一个阳性和一个阴性，模型给阳性打分更高的概率"。它只关心排序，不关心概率的绝对值准不准，所以是不平衡数据下相对稳健的指标（但也别过度迷信）。
library(pROC) r <- roc(df$y, phat, quiet = TRUE, direction = "<") auc(r) ci.auc(r) plot(r, print.auc = TRUE, main = "ROC curve")AUC的经验判读：0.5等于瞎猜，0.6~0.7一般，0.7~0.8尚可，0.8~0.9不错，0.9以上要么模型真强，要么有信息泄漏（比如把结局相关的字段误放进了特征里）。我遇到AUC突然飙到0.95以上时，第一反应不是高兴，是去查有没有泄漏。
要注意训练集上的AUC永远是乐观的，变量越多越乐观。所以交付前至少做一次交叉验证：
library(caret) set.seed(2024) ctrl <- trainControl(method = "cv", number = 10, classProbs = TRUE, summaryFunction = twoClassSummary, savePredictions = TRUE) fit_cv <- train(y ~ age + sex + bmi + smoke, data = df, method = "glm", family = "binomial", trControl = ctrl, metric = "ROC") fit_cv$results如果训练集AUC 0.85、交叉验证AUC 0.72，说明过拟合明显，该减变量或者加正则化。这个差距比单独看任何一个数都有信息量。
5.2 校准度与Hosmer-Lemeshow检验
AUC高不代表概率准。一个模型可以很好地把高低风险排开，但预测概率系统性偏高50%——这种情况下如果拿概率去做定价或者资源分配，会出大问题。
校准的检查分两步。先画校准曲线，按预测概率分十分位，比较每组的预测均值与实际事件率：
library(ResourceSelection) hoslem.test(df$y, phat, g = 10) # 或画校准图 library(rms) val.prob(phat, df$y)Hosmer-Lemeshow检验的原假设是"预测与观测吻合"，所以p值大才是好事，这一点跟大多数检验相反，容易搞混。而且它对分组数g很敏感，样本量大时几乎必然显著，小样本时又没效能。所以我只把它当参考，主要看校准曲线的形状。
如果发现校准不佳，有两种常见修正：一是重新校准（用Platt scaling或者等渗回归在验证集上拟合一个映射），二是检查是不是漏了重要变量或者函数形式不对。实践中，最常见的原因是连续变量的线性假设不成立，加个样条往往校准就改善了。这个观察我验证过很多次，值得优先排查。
5.3 内部验证：bootstrap与交叉验证
小样本项目里，我最推荐bootstrap做内部验证，因为它对样本量的浪费少：
library(rms) dd <- datadist(df); options(datadist = "dd") fit_lrm <- lrm(y ~ age + sex + bmi + smoke, data = df, x = TRUE, y = TRUE) validate(fit_lrm, method = "boot", B = 200)输出里会给出校正后的AUC、校准斜率等指标，还有" optimism"那一列，告诉你原始指标高估了多少。这个校正量在小样本下往往很可观，我在几十例的项目里见过AUC从0.82校正到0.71的情况。
交叉验证的好处是流程通用、容易解释，caret和tidymodels都能做。要注意的是，所有涉及数据的信息都必须在每一折内部重新计算，包括缺失值插补、标准化、变量筛选。一旦在划分折之前用了全量数据做标准化，就是在泄漏，AUC会虚高。这个坑极其常见，我自己也踩过。
# tidymodels 的写法，prep在每折内部执行 library(tidymodels) rec <- recipe(y ~ ., data = df) %>% step_normalize(all_numeric_predictors()) %>% step_impute_median(all_numeric_predictors())6. 踩坑记录与排查表
6.1 完全分离与准完全分离
这是logistic回归最经典的疑难杂症。现象是：summary里某个变量的系数特别大（比如5以上）、标准误更大（比如几十）、p值接近1，同时伴有fitted probabilities numerically 0 or 1 occurred的警告。原因通常是某个自变量能把结局完美区分开——比如某个二分类特征，阳性样本里有这个特征的100%是事件、0%是非事件。
这种情况下极大似然估计根本不存在（或者说趋于无穷），软件算出来的那个大数字没有意义。处理办法有三个：
第一，检查数据。先确认这不是编码错误或者泄漏。我遇到过"是否已取消"这种字段被误放进模型的，它当然能完美预测"是否流失"，因为它是事后信息。
第二，用Firth惩罚似然。这是小样本和分离问题的标准解法，通过给似然加一个惩罚项来消除一阶偏差：
library(logistf) fit_firth <- logistf(y ~ age + sex + bmi, data = df) summary(fit_firth)第三，合并类别或者干脆删掉这个变量。如果一个变量分离得太厉害，它对模型的贡献是"确定性规则"，而这种规则往往泛化性很差。
准分离（quasi-separation）没那么极端，但也会让系数膨胀。判断方法是看这个变量分出来的各组事件数，如果某一组只有几个事件甚至0个，基本可以确定。
注意：分离问题在小样本、稀疏类别、强预测变量三个条件凑齐时最容易出现。建模前的交叉表检查能提前发现：对每个分类变量做
table(df$var, df$y)，看有没有单元格是0。
6.2 收敛失败与警告信息
algorithm did not converge这个警告我一年能碰上七八次，原因无非几类：
数值尺度问题最常见。一个变量取值上百万，另一个0~1，Fisher scoring迭代时梯度差别太大。解决办法就是前面说的标准化。
初始值问题。默认从0开始，某些情况下换一个起点能收敛：
fit <- glm(y ~ ., data = df, family = binomial, start = rep(0, ncol(model.matrix(y ~ ., df))))共线性太严重也会导致不收敛。这时候看vif，把VIF最大的那个拿掉试试。
如果真的怎么调都不收敛，考虑换算法：glm2包的glm2()在数值稳定性上比glm好一些，接口基本兼容。或者干脆用正则化方法，glmnet天然能处理这类问题。
library(glm2) fit2 <- glm2(y ~ ., data = df, family = binomial)另外提醒一句：警告信息不要用suppressWarnings()一盖了之。那等于把仪表盘上的故障灯贴了张纸。
6.3 常见问题速查表
下面这张表是我自己攒的，遇到问题先扫一眼，能省不少时间：
| 现象 | 可能原因 | 排查动作 | 处理方案 |
|---|---|---|---|
| OR值极大且区间极宽 | 完全/准分离 | table(var, y)看零单元格 | Firth惩罚似然、合并类别、删变量 |
| 多个变量p值都很大 | 共线性 | vif(fit) | 合并变量、正则化、删冗余 |
| 不收敛警告 | 尺度差异、共线、分离 | 标准化后再试 | 标准化、glm2、加正则 |
| 训练AUC高、CV AUC低 | 过拟合、变量过多 | 对比两者差值 | 减变量、加惩罚、增加样本 |
| 校准曲线偏离对角线 | 线性假设不成立 | 连续变量分箱看趋势 | 加样条、分段、变换 |
| 系数符号与常识相反 | 共线性、抑制效应、编码 | 单独建单变量模型对比 | 检查编码、看相关矩阵 |
| 预测概率全是0或1 | 分离、尺度未处理 | 看phat分布 | 同上，优先查分离 |
| 因子变量的OR不知对应哪组 | 参照水平搞错 | levels(df$var) | relevel()显式指定参照 |
再说两个容易被忽略的细节。第一个是因子变量的参照水平：relevel一行就能改，但如果不显式指定，R默认按字母序，报告里每个OR对应哪一组就容易说不清。第二个是样本内预测和样本外预测的差距：如果新数据的预测分布和训练集差别很大（比如训练集里概率大多在0.1~0.3，新数据里大多在0.6以上），先别怀疑模型，去检查新数据的变量分布和编码是不是和训练时一致——这类"数据管道错位"造成的问题，比模型本身的问题多得多。
最后分享一个我自己的习惯：每跑完一个模型，我会用saveRDS(fit, "model.rds")把模型对象存下来，同时把建模用的数据版本、变量列表、R包版本记在一个小文本里。听起来很琐碎，但当下个月有人问"这个模型当时用的哪些变量"时，你能在三秒内答上来，这种踏实感是省不出来的。下一篇我会接着聊变量筛选的具体策略、非线性关系的处理，以及怎么把模型结果包装成业务方愿意看的形式。
