The V Lab
适用对象公共卫生、流行病学及健康科学学习者
学习时长约 120–180 分钟
先修要求概率、置信区间与基础 R

关于本教程的数据 全部记录均由固定随机种子模拟生成，不包含真实个人健康信息。模拟机制仅用于教学；代码中的关联不应被理解为真实世界中的因果效应或临床效应。

如何使用本教程

本教程以“是否发生呼吸道感染”这一二元结局为主线。建议先依次阅读概率、优势和 logit，再运行模型代码；随后重点练习把模型结果翻译回概率尺度。所有示例仅依赖 R 自带的 stats、graphics 和 knitr，可以在干净的 R 会话中 Knit。

学习目标

完成本教程后，你应该能够：

  1. 判断研究问题是否适合逻辑回归；
  2. 区分概率（probability）、优势（odds）、对数优势（log-odds）和优势比（odds ratio, OR）；
  3. 使用 glm(..., family = binomial()) 拟合并解释粗模型和调整模型；
  4. 将模型系数转换为 OR、置信区间和预测概率；
  5. 处理分类变量、非线性关系和交互作用；
  6. 区分校准、判别和阈值分类，并在测试集上评价模型；
  7. 检查残差、杠杆值、影响点、稀疏数据和分离问题；
  8. 以透明、不过度因果化的方式报告结果。

1 为什么需要逻辑回归

1.1 二元结局不能直接套用普通线性回归

逻辑回归适用于每个观测的结局只能取两个互斥状态的情形，例如：

  • 感染 / 未感染；
  • 死亡 / 存活；
  • 筛查阳性 / 阴性；
  • 接受服务 / 未接受服务。

若把 0/1 结局直接放进普通线性回归，预测值可能小于 0 或大于 1，误差方差也会随均值变化。逻辑回归通过链接函数把任意实数映射到 0 与 1 之间，从而对事件概率建模。

先定义“事件” 模型开始前必须写明哪个值代表事件。本教程中 infection_num = 1 表示“发生感染”。因子的水平顺序也显式设为“否”“是”，但核心模型使用数值型 0/1 变量，避免事件方向含糊。

1.2 概率、优势与 logit

若事件概率为 (p)，则：

odds=p1−p,logit(p)=log⁡(p1−p). \text{odds}=\frac{p}{1-p}, \qquad \text{logit}(p)=\log\left(\frac{p}{1-p}\right).

概率是“事件数 / 总人数”；优势是“事件数 / 非事件数”。二者数值通常不同。

probability_examples <- c(0.05, 0.20, 0.50, 0.80)
probability_table <- data.frame(
  概率 = probability_examples,
  优势 = probability_examples / (1 - probability_examples),
  对数优势 = qlogis(probability_examples)
)
knitr::kable(
  probability_table,
  digits = 3,
  caption = "概率、优势和对数优势的对应关系"
)
概率、优势和对数优势的对应关系
概率 优势 对数优势
0.05 0.053 -2.94
0.20 0.250 -1.39
0.50 1.000 0.00
0.80 4.000 1.39

当 (p=0.20) 时，优势为 (0.20/0.80=0.25)，即大约 1:4；不能把 0.25 说成 25% 的事件概率。反向转换使用：

p=odds1+odds=exp⁡(η)1+exp⁡(η). p=\frac{\text{odds}}{1+\text{odds}} =\frac{\exp(\eta)}{1+\exp(\eta)}.

1.3 模型形式

含 (k) 个预测变量的逻辑回归写为：

log⁡(pi1−pi)=β0+β1Xi1+⋯+βkXik. \log\left(\frac{p_i}{1-p_i}\right) =\beta_0+\beta_1X_{i1}+\cdots+\beta_kX_{ik}.

模型假设预测变量与对数优势呈指定的函数关系，而不是假设概率与预测变量呈直线关系。系数通过最大似然法估计：在候选参数中，寻找使已观察到的 0/1 结局最可能出现的一组参数。

2 认识数据并提出问题

2.1 研究问题与变量角色

示例问题是：

在这份模拟数据中，接种状态与一年内呼吸道感染是否相关？在控制年龄、BMI、吸烟和居住地区后，这种关联如何？模型对新观测的概率预测表现如何？

这里的结局为感染，主要解释变量为接种状态，其他变量可能用于减少混杂或改善预测。变量是否应调整，应由研究问题、时间顺序和领域知识决定，而不是由单变量 p 值筛选。

模拟完整数据已在任何结局探索之前分层保留 30% 作为测试集。下面的描述、探索、函数形式选择和模型拟合全部只使用 train_data；test_data 要到最终评价一节才首次汇总。

str(train_data)
## 'data.frame':    839 obs. of  10 variables:
##  $ participant_id: chr  "P0001" "P0002" "P0003" "P0004" ...
##  $ infection_num : int  0 0 0 0 0 0 0 0 1 0 ...
##  $ infection     : Factor w/ 2 levels "否","是": 1 1 1 1 1 1 1 1 2 1 ...
##  $ age           : num  26 32 61 33 40 32 59 54 41 41 ...
##  $ age10         : num  -2.4 -1.8 1.1 -1.7 -1 -1.8 0.9 0.4 -0.9 -0.9 ...
##  $ bmi           : num  31.5 31.1 28 30.8 30.3 31.9 23.8 27.3 28.8 23 ...
##  $ bmi5          : num  1.3 1.22 0.6 1.16 1.06 1.38 -0.24 0.46 0.76 -0.4 ...
##  $ smoking       : Factor w/ 2 levels "否","是": 2 1 1 1 1 1 2 1 1 1 ...
##  $ vaccinated    : Factor w/ 2 levels "否","是": 2 1 2 1 1 1 2 2 1 2 ...
##  $ area          : Factor w/ 2 levels "城市","农村": 2 1 1 1 1 2 1 2 1 1 ...
data_summary <- data.frame(
  样本量 = nrow(train_data),
  感染人数 = sum(train_data$infection_num),
  感染比例 = mean(train_data$infection_num),
  平均年龄 = mean(train_data$age),
  平均BMI = mean(train_data$bmi)
)
knitr::kable(data_summary, digits = 3, caption = "开发/训练数据概览")
开发/训练数据概览
样本量 感染人数 感染比例 平均年龄 平均BMI
839 169 0.201 48.6 26.4

2.2 先看绝对人数和风险

vaccination_table <- with(
  train_data,
  table(接种 = vaccinated, 感染 = infection)
)
vaccination_table
##     感染
## 接种  否  是
##   否 310  98
##   是 360  71
infection_by_vaccination <- aggregate(
  infection_num ~ vaccinated,
  data = train_data,
  FUN = function(x) c(n = length(x), events = sum(x), risk = mean(x))
)
infection_by_vaccination <- data.frame(
  接种 = infection_by_vaccination$vaccinated,
  样本量 = infection_by_vaccination$infection_num[, "n"],
  感染人数 = infection_by_vaccination$infection_num[, "events"],
  感染比例 = infection_by_vaccination$infection_num[, "risk"]
)
knitr::kable(
  infection_by_vaccination,
  digits = 3,
  caption = "按接种状态汇总的感染人数与比例"
)
按接种状态汇总的感染人数与比例
接种 样本量 感染人数 感染比例
否 408 98 0.240
是 431 71 0.165 
observed_risk <- tapply(
  train_data$infection_num,
  train_data$vaccinated,
  mean
)
group_n <- table(train_data$vaccinated)
observed_se <- sqrt(observed_risk * (1 - observed_risk) / group_n)

plot(
  seq_along(observed_risk), observed_risk,
  pch = c(16, 17), cex = 1.3,
  col = c(palette_lr["vermillion"], palette_lr["blue"]),
  xaxt = "n", ylim = c(0, max(observed_risk + 2 * observed_se) * 1.10),
  xlab = "接种状态", ylab = "观察感染比例",
  main = "先查看绝对风险"
)
axis(1, at = seq_along(observed_risk), labels = names(observed_risk))
arrows(
  seq_along(observed_risk), observed_risk - 1.96 * observed_se,
  seq_along(observed_risk), observed_risk + 1.96 * observed_se,
  angle = 90, code = 3, length = 0.06,
  col = c(palette_lr["vermillion"], palette_lr["blue"])
)
点图比较未接种组和接种组的观察感染比例，并带有垂直误差线。

不同接种组的观察感染比例；误差线为近似 95% 置信区间。

描述性比较很重要，但它尚未控制各组构成差异，也不能自动解释为接种造成的效果。

3 拟合第一个逻辑回归模型

3.1 粗模型

crude_model <- glm(
  infection_num ~ vaccinated,
  data = train_data,
  family = binomial(link = "logit")
)
summary(crude_model)
## 
## Call:
## glm(formula = infection_num ~ vaccinated, family = binomial(link = "logit"), 
##     data = train_data)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -1.152      0.116   -9.94   <2e-16 ***
## vaccinated是   -0.472      0.174   -2.71   0.0067 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 842.99  on 838  degrees of freedom
## Residual deviance: 835.56  on 837  degrees of freedom
## AIC: 839.6
## 
## Number of Fisher Scoring iterations: 4

glm() 默认输出 log-odds 尺度上的系数。对接种系数取指数，得到接种者相对于未接种者的粗 OR。

crude_beta <- coef(crude_model)["vaccinated是"]
crude_se <- summary(crude_model)$coefficients["vaccinated是", "Std. Error"]
crude_or <- exp(crude_beta)
crude_ci <- exp(crude_beta + qnorm(c(0.025, 0.975)) * crude_se)

crude_result <- data.frame(
  对比 = "已接种 vs 未接种",
  OR = crude_or,
  置信区间下限 = crude_ci[1],
  置信区间上限 = crude_ci[2]
)
knitr::kable(crude_result, digits = 3, caption = "接种状态的粗优势比")
接种状态的粗优势比
对比 OR 置信区间下限 置信区间上限
vaccinated是 已接种 vs 未接种 0.624 0.444 0.877

训练数据中，接种者相对于未接种者的感染优势比为 0.62（Wald 95% 置信区间：0.44 至 0.88）。这描述的是优势之比，不是概率之比，也不是风险降低百分比。

3.2 调整模型与有意义的单位

年龄按 10 年、BMI 按 5 kg/m² 表示，并分别以 50 岁和 BMI 25 为中心。这样截距有现实含义，连续变量 OR 的单位也更容易沟通。

adjusted_model <- glm(
  infection_num ~ age10 + bmi5 + smoking + vaccinated + area,
  data = train_data,
  family = binomial()
)

model_coefficient_table <- function(model, labels = NULL) {
  coefficient_matrix <- summary(model)$coefficients
  output <- data.frame(
    项 = rownames(coefficient_matrix),
    系数 = coefficient_matrix[, "Estimate"],
    标准误 = coefficient_matrix[, "Std. Error"],
    OR = exp(coefficient_matrix[, "Estimate"]),
    OR下限 = exp(
      coefficient_matrix[, "Estimate"] -
        qnorm(0.975) * coefficient_matrix[, "Std. Error"]
    ),
    OR上限 = exp(
      coefficient_matrix[, "Estimate"] +
        qnorm(0.975) * coefficient_matrix[, "Std. Error"]
    ),
    p值 = coefficient_matrix[, "Pr(>|z|)"],
    row.names = NULL
  )
  if (!is.null(labels)) {
    replacement <- unname(labels[output$项])
    output$项[!is.na(replacement)] <- replacement[!is.na(replacement)]
  }
  output
}

adjusted_labels <- c(
  "(Intercept)" = "截距：50岁、BMI 25、非吸烟、未接种、城市",
  age10 = "年龄（每增加10年）",
  bmi5 = "BMI（每增加5 kg/m²）",
  "smoking是" = "吸烟：是 vs 否",
  "vaccinated是" = "接种：是 vs 否",
  "area农村" = "地区：农村 vs 城市"
)

adjusted_results <- model_coefficient_table(adjusted_model, adjusted_labels)
knitr::kable(
  adjusted_results,
  digits = 3,
  caption = "感染结局的多变量逻辑回归（Wald 置信区间）"
)
感染结局的多变量逻辑回归（Wald 置信区间）
项 系数 标准误 OR OR下限 OR上限 p值
截距：50岁、BMI 25、非吸烟、未接种、城市 -1.600 0.158 0.202 0.148 0.275 0.000
年龄（每增加10年） 0.424 0.065 1.527 1.344 1.735 0.000
BMI（每增加5 kg/m²） 0.405 0.105 1.499 1.219 1.843 0.000
吸烟：是 vs 否 0.748 0.208 2.112 1.406 3.173 0.000
接种：是 vs 否 -0.652 0.188 0.521 0.361 0.753 0.001
地区：农村 vs 城市 0.468 0.190 1.596 1.100 2.315 0.014

对于接种系数，解释必须包含“在模型中其他变量相同”的条件；对于年龄和 BMI，必须写明单位。

在年龄、BMI、吸烟状态和地区相同的观测之间，接种者相对于未接种者的估计感染优势比为 0.52。年龄每增加 10 年的调整后 OR 为 1.53。这些是条件关联；模拟的观测性分析本身不能证明因果关系。

3.2.1 截距应该怎样理解

截距是所有数值变量等于 0、所有分类变量位于参照水平时的 log-odds。由于年龄和 BMI 已中心化，这里对应 50 岁、BMI 25、非吸烟、未接种且居住城市的人。截距取指数是该参照画像的基线优势，不是某个比较的 OR。

3.3 OR 为什么不是风险比

若未暴露组风险为 (p_0)，OR 为 ( heta)，则暴露组对应风险为：

p1=θp01−p0+θp0. p_1=\frac{\theta p_0}{1-p_0+\theta p_0}. 

baseline_risks <- c(0.05, 0.20, 0.50)
example_or <- 2
risk_from_or <- example_or * baseline_risks /
  (1 - baseline_risks + example_or * baseline_risks)

or_risk_table <- data.frame(
  基线风险 = baseline_risks,
  OR = example_or,
  对应风险 = risk_from_or,
  风险比 = risk_from_or / baseline_risks,
  风险差 = risk_from_or - baseline_risks
)
knitr::kable(
  or_risk_table,
  digits = 3,
  caption = "同一个 OR 在不同基线风险下对应不同风险比和风险差"
)
同一个 OR 在不同基线风险下对应不同风险比和风险差
基线风险 OR 对应风险 风险比 风险差
0.05 2 0.095 1.91 0.045
0.20 2 0.333 1.67 0.133
0.50 2 0.667 1.33 0.167

只有在结局较罕见时，OR 与风险比才可能数值接近。即使如此，报告时仍应使用正确名称。

4 从系数回到预测概率

4.1 为具体画像预测

profile_data <- expand.grid(
  age10 = c(0, 2),
  bmi5 = 0,
  smoking = factor("否", levels = levels(train_data$smoking)),
  vaccinated = factor(c("否", "是"), levels = levels(train_data$vaccinated)),
  area = factor("城市", levels = levels(train_data$area))
)
profile_data$年龄 <- profile_data$age10 * 10 + 50
profile_data$BMI <- profile_data$bmi5 * 5 + 25

link_prediction <- predict(
  adjusted_model,
  newdata = profile_data,
  type = "link",
  se.fit = TRUE
)
profile_data$预测概率 <- plogis(link_prediction$fit)
profile_data$概率下限 <- plogis(
  link_prediction$fit - qnorm(0.975) * link_prediction$se.fit
)
profile_data$概率上限 <- plogis(
  link_prediction$fit + qnorm(0.975) * link_prediction$se.fit
)

knitr::kable(
  profile_data[, c(
    "年龄", "BMI", "smoking", "vaccinated", "area",
    "预测概率", "概率下限", "概率上限"
  )],
  digits = 3,
  col.names = c(
    "年龄", "BMI", "吸烟", "接种", "地区",
    "预测概率", "均值概率95%下限", "均值概率95%上限"
  ),
  caption = "具体协变量画像的模型预测概率"
)
具体协变量画像的模型预测概率
年龄 BMI 吸烟 接种 地区 预测概率 均值概率95%下限 均值概率95%上限
50 25 否 否 城市 0.168 0.129 0.216
70 25 否 否 城市 0.320 0.242 0.410
50 25 否 是 城市 0.095 0.070 0.129
70 25 否 是 城市 0.197 0.144 0.263

应先在 link 尺度构造置信区间，再用 plogis() 转回概率尺度。这里的区间反映该画像平均概率估计的不确定性，不是对某个个体最终 0/1 结局的概率区间。

4.2 标准化预测与平均风险差

OR 往往不直观。可以把训练数据中的每个人复制两次，一次设为未接种，一次设为已接种，再分别平均模型预测概率。这称为模型标准化或预测边际化。

standardized_no <- train_data
standardized_yes <- train_data
standardized_no$vaccinated <- factor(
  "否", levels = levels(train_data$vaccinated)
)
standardized_yes$vaccinated <- factor(
  "是", levels = levels(train_data$vaccinated)
)

standardized_risk_no <- mean(predict(
  adjusted_model, newdata = standardized_no, type = "response"
))
standardized_risk_yes <- mean(predict(
  adjusted_model, newdata = standardized_yes, type = "response"
))

standardized_results <- data.frame(
  情景 = c("所有人设为未接种", "所有人设为已接种", "已接种 - 未接种"),
  估计 = c(
    standardized_risk_no,
    standardized_risk_yes,
    standardized_risk_yes - standardized_risk_no
  )
)
knitr::kable(
  standardized_results,
  digits = 3,
  caption = "基于调整模型的标准化概率和平均风险差"
)
基于调整模型的标准化概率和平均风险差
情景 估计
所有人设为未接种 0.250
所有人设为已接种 0.158
已接种 - 未接种 -0.093

标准化并不会自动产生因果效应 把预测变量人为设为两个水平是一种计算方式。若要把差值解释为因果效应，还需要一致性、可交换性、正值性、正确模型形式和可靠测量等额外假设。

5 分类变量、非线性与交互作用

5.1 分类变量与参照水平

R 对含两个水平的因子建立一个虚拟变量。系数比较非参照水平和参照水平。先显式检查水平：

lapply(
  train_data[c("smoking", "vaccinated", "area", "infection")],
  levels
)
## $smoking
## [1] "否" "是"
## 
## $vaccinated
## [1] "否" "是"
## 
## $area
## [1] "城市" "农村"
## 
## $infection
## [1] "否" "是"
model.matrix(
  ~ smoking + vaccinated + area,
  data = train_data[1:6, ]
)
##    (Intercept) smoking是 vaccinated是 area农村
## 1            1         1            1        1
## 2            1         0            0        0
## 3            1         0            1        0
## 4            1         0            0        0
## 5            1         0            0        0
## 10           1         0            0        1
## attr(,"assign")
## [1] 0 1 2 3
## attr(,"contrasts")
## attr(,"contrasts")$smoking
## [1] "contr.treatment"
## 
## attr(,"contrasts")$vaccinated
## [1] "contr.treatment"
## 
## attr(,"contrasts")$area
## [1] "contr.treatment"

如需改变参照组，可用 relevel()。改变参照组会改变系数写法，却不会改变每个观测的拟合概率。

5.2 连续预测变量在 logit 尺度上的线性

逻辑回归不要求年龄本身服从正态分布，但基础模型假设年龄与 logit 之间为直线关系。直接把未调整的年龄组感染率与某个固定协变量画像的条件预测作比较，会混合不同目标量；这里改用成分加残差图（component-plus-residual plot），在调整模型的同一条件 logit 尺度上检查年龄项。

age_component <-
  unname(coef(adjusted_model)["age10"]) * train_data$age10
age_partial_residual <-
  age_component + residuals(adjusted_model, type = "working")

age_group <- cut(
  train_data$age,
  breaks = quantile(train_data$age, probs = seq(0, 1, 0.1)),
  include.lowest = TRUE,
  ordered_result = TRUE
)
age_partial_data <- data.frame(
  age = train_data$age,
  partial_residual = age_partial_residual,
  age_group = age_group
)
age_partial_summary <- aggregate(
  cbind(age, partial_residual) ~ age_group,
  data = age_partial_data,
  FUN = mean
)

age_sequence <- seq(min(train_data$age), max(train_data$age), length.out = 200)
linear_age_component <-
  unname(coef(adjusted_model)["age10"]) * ((age_sequence - 50) / 10)
age_smooth <- lowess(
  train_data$age,
  age_partial_residual,
  f = 2 / 3
)
plot(
  train_data$age, age_partial_residual,
  pch = 16, cex = 0.45,
  col = rgb(107/255, 114/255, 128/255, 0.28),
  xlab = "年龄（岁）", ylab = "年龄项成分 + 工作残差",
  main = "在调整后的条件 logit 尺度检查年龄函数形式"
)
points(
  age_partial_summary$age,
  age_partial_summary$partial_residual,
  pch = 16, cex = 1.15, col = palette_lr["blue"]
)
lines(
  age_smooth$x, age_smooth$y,
  lwd = 3, lty = 2, col = palette_lr["teal"]
)
lines(
  age_sequence, linear_age_component,
  lwd = 3, col = palette_lr["vermillion"]
)
legend(
  "topleft",
  legend = c("年龄十分位均值", "LOWESS 平滑", "模型指定的线性年龄项"),
  pch = c(16, NA, NA), lty = c(NA, 2, 1), lwd = c(NA, 3, 3),
  col = c(
    palette_lr["blue"],
    palette_lr["teal"],
    palette_lr["vermillion"]
  ),
  bty = "n"
)
散点和平滑图在调整后的条件 logit 尺度上比较年龄部分残差与模型指定的线性年龄项。

调整模型中年龄的成分加残差图。分组均值和平滑线若系统偏离模型直线，提示应重新考虑函数形式。

工作残差在预测概率接近 0 或 1 时可能很大，因此图形只是一种诊断线索。若平滑线持续偏离模型直线，应结合领域知识比较预先指定的平方项、分段线性项或样条，并用重采样评价是否真正改善预测。

若领域知识或图形提示弯曲关系，可以预先指定平方项、分段线性项或限制性立方样条。下面比较线性年龄项和平方年龄项；这只是函数形式演示，不应以不断试模来追逐较小 p 值。

quadratic_model <- glm(
  infection_num ~ age10 + I(age10^2) + bmi5 + smoking + vaccinated + area,
  data = train_data,
  family = binomial()
)

anova(adjusted_model, quadratic_model, test = "Chisq")
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
833 748 NA NA NA
832 748 1 0.004 0.952 
AIC(adjusted_model, quadratic_model)
df AIC
adjusted_model 6 760
quadratic_model 7 762

5.3 交互作用：效应取决于情境

若接种关联可能因地区而异，可以加入接种与地区的乘积项：

interaction_model <- glm(
  infection_num ~ age10 + bmi5 + smoking + vaccinated * area,
  data = train_data,
  family = binomial()
)

interaction_coefficients <- coef(interaction_model)
interaction_vcov <- vcov(interaction_model)

# 城市中，接种 OR 只涉及 vaccinated是 系数。
urban_contrast <- c(
  "vaccinated是" = 1,
  "vaccinated是:area农村" = 0
)

# 农村中，接种 log(OR) 是主效应与交互项之和。
rural_contrast <- c(
  "vaccinated是" = 1,
  "vaccinated是:area农村" = 1
)

contrast_or <- function(model, contrast) {
  all_contrast <- setNames(rep(0, length(coef(model))), names(coef(model)))
  all_contrast[names(contrast)] <- contrast
  estimate <- sum(all_contrast * coef(model))
  standard_error <- sqrt(
    as.numeric(t(all_contrast) %*% vcov(model) %*% all_contrast)
  )
  c(
    OR = exp(estimate),
    lower = exp(estimate - qnorm(0.975) * standard_error),
    upper = exp(estimate + qnorm(0.975) * standard_error)
  )
}

area_specific_or <- rbind(
  城市 = contrast_or(interaction_model, urban_contrast),
  农村 = contrast_or(interaction_model, rural_contrast)
)
knitr::kable(
  data.frame(地区 = rownames(area_specific_or), area_specific_or, row.names = NULL),
  digits = 3,
  col.names = c("地区", "接种OR", "95%下限", "95%上限"),
  caption = "交互模型中按地区计算的接种优势比"
)
交互模型中按地区计算的接种优势比
地区 接种OR 95%下限 95%上限
城市 0.406 0.254 0.65
农村 0.780 0.433 1.40

含交互项时，vaccinated是 只表示参照地区“城市”中的接种比较；农村中的接种比较必须再加上交互项。不要孤立解读所谓“主效应”。

6 模型拟合、比较与诊断

6.1 似然、偏差与伪 R²

较低的残差偏差和 AIC 可以帮助比较在同一数据、同一结局上拟合的候选模型。嵌套模型可用似然比检验。它们不能替代外部验证或科学合理性。

null_model <- glm(
  infection_num ~ 1,
  data = train_data,
  family = binomial()
)

mcfadden_r2 <- 1 - as.numeric(logLik(adjusted_model) / logLik(null_model))
fit_statistics <- data.frame(
  指标 = c("空模型偏差", "调整模型偏差", "AIC", "McFadden伪R²"),
  数值 = c(
    deviance(null_model),
    deviance(adjusted_model),
    AIC(adjusted_model),
    mcfadden_r2
  )
)
knitr::kable(fit_statistics, digits = 3, caption = "训练数据的拟合指标")
训练数据的拟合指标
指标 数值
空模型偏差 842.992
调整模型偏差 748.125
AIC 760.125
McFadden伪R² 0.113 
anova(null_model, adjusted_model, test = "Chisq")
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
838 843 NA NA NA
833 748 5 94.9 0

McFadden 伪 R² 不是线性回归中“解释方差比例”的 (R^2)，数值也不能直接比较。报告时应说明定义。

6.2 逻辑回归不要求残差正态或同方差

逻辑回归的响应是伯努利变量，其方差天然为 (p(1-p))。因此不应把线性回归的残差正态性和同方差性要求照搬过来。应检查函数形式、异常残差、杠杆值、影响点、相关观测、过度离散、分离和模型设定。

diagnostic_data <- data.frame(
  fitted_probability = fitted(adjusted_model),
  deviance_residual = residuals(adjusted_model, type = "deviance"),
  pearson_residual = residuals(adjusted_model, type = "pearson"),
  leverage = hatvalues(adjusted_model),
  cooks_distance = cooks.distance(adjusted_model)
)

diagnostic_summary <- data.frame(
  指标 = c(
    "|偏差残差| > 2",
    "杠杆值 > 2p/n",
    "Cook距离 > 4/n"
  ),
  数量 = c(
    sum(abs(diagnostic_data$deviance_residual) > 2),
    sum(
      diagnostic_data$leverage >
        2 * length(coef(adjusted_model)) / nrow(train_data)
    ),
    sum(diagnostic_data$cooks_distance > 4 / nrow(train_data))
  )
)
knitr::kable(
  diagnostic_summary,
  caption = "用于定位需进一步核查观测的经验阈值"
)
用于定位需进一步核查观测的经验阈值
指标 数量
|偏差残差| > 2 25
杠杆值 > 2p/n 72
Cook距离 > 4/n 59 
old_par <- par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))

plot(
  diagnostic_data$fitted_probability,
  diagnostic_data$deviance_residual,
  pch = 16, cex = 0.55,
  col = rgb(0, 114/255, 178/255, 0.40),
  xlab = "拟合概率", ylab = "偏差残差",
  main = "残差"
)
abline(h = c(-2, 0, 2), lty = c(2, 1, 2), col = palette_lr["grey"])

plot(
  diagnostic_data$leverage,
  type = "h", col = palette_lr["teal"],
  xlab = "训练集观测序号", ylab = "杠杆值",
  main = "杠杆值"
)
abline(
  h = 2 * length(coef(adjusted_model)) / nrow(train_data),
  lty = 2, col = palette_lr["vermillion"]
)

plot(
  diagnostic_data$cooks_distance,
  type = "h", col = palette_lr["orange"],
  xlab = "训练集观测序号", ylab = "Cook 距离",
  main = "影响度"
)
abline(
  h = 4 / nrow(train_data),
  lty = 2, col = palette_lr["vermillion"]
)
三联图分别显示拟合概率与偏差残差、观测序号与杠杆值、观测序号与 Cook 距离。

逻辑回归的偏差残差、杠杆值和 Cook 距离诊断。经验阈值用于筛查，而不是自动删除规则。 

par(old_par)

高残差表示结局与模型预测不一致；高杠杆值表示协变量组合少见；高 Cook 距离表示删除该观测可能较明显地改变拟合。应回查数据质量、研究过程和敏感性分析，而不是看到阈值就机械删除。

6.3 稀疏单元、完全分离与收敛

若某个预测变量水平中所有人都发生事件或都不发生事件，最大似然估计可能趋向正负无穷，这称为完全分离。常见信号包括：

  • 交叉表有 0 单元；
  • 系数绝对值和标准误极大；
  • glm() 给出“拟合概率为 0 或 1”或未收敛警告；
  • 不同软件给出极不稳定结果。
with(train_data, table(infection, smoking, area))
## , , area = 城市
## 
##          smoking
## infection  否  是
##        否 407  76
##        是  69  30
## 
## , , area = 农村
## 
##          smoking
## infection  否  是
##        否 151  36
##        是  45  25
c(
  converged = adjusted_model$converged,
  iterations = adjusted_model$iter,
  events = sum(train_data$infection_num),
  non_events = sum(train_data$infection_num == 0),
  parameters = length(coef(adjusted_model))
)
##  converged iterations     events non_events parameters 
##          1          5        169        670          6

可考虑减少不必要参数、合并有科学依据的类别、增加有信息的样本，或使用惩罚/偏倚校正方法。固定的“每个参数 10 个事件”只是过度简化的经验规则；事件数、非事件数、效应大小、预测变量分布、缺失和验证方式都影响可靠性。

6.4 过度离散与相关观测

对独立的个体级二元数据，理论上伯努利方差已经由均值决定。若数据是分组计数、重复测量、家庭或社区聚类，简单逻辑回归的独立性和方差设定可能不成立。此时可考虑准二项、广义估计方程或混合效应模型，并按抽样设计估计不确定性。

7 预测性能：校准与判别

7.1 必须在未用于拟合的数据上评估

训练集用于估计系数，测试集只用于最终评价。这里的 70/30 分割按事件状态分层，因此两个集合都包含事件和非事件。

split_summary <- data.frame(
  数据集 = c("训练集", "测试集"),
  样本量 = c(nrow(train_data), nrow(test_data)),
  事件数 = c(sum(train_data$infection_num), sum(test_data$infection_num)),
  事件比例 = c(mean(train_data$infection_num), mean(test_data$infection_num))
)
knitr::kable(split_summary, digits = 3, caption = "训练集与测试集构成")
训练集与测试集构成
数据集 样本量 事件数 事件比例
训练集 839 169 0.201
测试集 361 73 0.202 
test_probability <- predict(
  adjusted_model,
  newdata = test_data,
  type = "response"
)

7.2 校准：预测概率是否与观察比例一致

Brier 分数是预测概率与 0/1 结局之间均方误差，越小越好；它同时受到校准和判别影响，并依赖结局基线率。

test_brier <- mean((test_data$infection_num - test_probability)^2)
test_log_loss <- -mean(
  test_data$infection_num * log(clamp_probability(test_probability)) +
    (1 - test_data$infection_num) *
      log(1 - clamp_probability(test_probability))
)

calibration_group <- cut(
  test_probability,
  breaks = unique(quantile(test_probability, probs = seq(0, 1, 0.1))),
  include.lowest = TRUE,
  ordered_result = TRUE
)
calibration_data <- aggregate(
  cbind(predicted = test_probability, observed = test_data$infection_num) ~
    calibration_group,
  FUN = mean
)

calibration_metrics <- data.frame(
  指标 = c("测试集 Brier 分数", "测试集 log loss"),
  数值 = c(test_brier, test_log_loss)
)
knitr::kable(calibration_metrics, digits = 3, caption = "测试集概率误差")
测试集概率误差
指标 数值
测试集 Brier 分数 0.156
测试集 log loss 0.486 
plot(
  calibration_data$predicted,
  calibration_data$observed,
  pch = 16, cex = 1.15, col = palette_lr["blue"],
  xlim = c(0, 1), ylim = c(0, 1), asp = 1,
  xlab = "组内平均预测概率", ylab = "组内观察事件比例",
  main = "测试集校准"
)
abline(0, 1, lty = 2, lwd = 2, col = palette_lr["grey"])
校准散点图比较各组平均预测概率和观察事件比例，并有一条理想的对角线。

测试集按预测概率十分位分组的校准图。理想状态接近 45 度线。

分组校准图会受分组方式和样本量影响，不应把某个校准检验的 p 值当作“模型通过”的证明。还应考察总体校准、关键概率范围和重要亚组。

7.3 判别：事件者是否得到更高分数

ROC 曲线展示所有阈值下灵敏度与 1−特异度的权衡。AUC 可解释为随机抽取一名事件者和一名非事件者时，模型给事件者更高分数的概率（并列计一半）。

auc_rank <- function(y, probability) {
  n_event <- sum(y == 1)
  n_nonevent <- sum(y == 0)
  rank_sum <- sum(rank(probability, ties.method = "average")[y == 1])
  (rank_sum - n_event * (n_event + 1) / 2) /
    (n_event * n_nonevent)
}

roc_coordinates <- function(y, probability) {
  thresholds <- c(Inf, sort(unique(probability), decreasing = TRUE), -Inf)
  true_positive_rate <- vapply(thresholds, function(threshold) {
    predicted <- probability >= threshold
    sum(predicted & y == 1) / sum(y == 1)
  }, numeric(1))
  false_positive_rate <- vapply(thresholds, function(threshold) {
    predicted <- probability >= threshold
    sum(predicted & y == 0) / sum(y == 0)
  }, numeric(1))
  data.frame(
    threshold = thresholds,
    false_positive_rate = false_positive_rate,
    true_positive_rate = true_positive_rate
  )
}

test_auc <- auc_rank(test_data$infection_num, test_probability)
test_roc <- roc_coordinates(test_data$infection_num, test_probability)
plot(
  test_roc$false_positive_rate,
  test_roc$true_positive_rate,
  type = "l", lwd = 3, col = palette_lr["blue"],
  xlim = c(0, 1), ylim = c(0, 1), asp = 1,
  xlab = "假阳性率（1 - 特异度）",
  ylab = "真阳性率（灵敏度）",
  main = paste0("测试集 ROC：AUC = ", round(test_auc, 3))
)
abline(0, 1, lty = 2, col = palette_lr["grey"])
ROC 曲线显示假阳性率与真阳性率之间的权衡，并以对角线表示随机排序。

调整模型在测试集上的 ROC 曲线。

AUC 为 0.648。AUC 衡量排序而不直接评价概率是否准确；一个 AUC 较高的模型仍可能系统性高估风险。

7.4 阈值、混淆矩阵和决策后果

把概率转成 0/1 分类必须选择阈值。0.5 并非天然最优；阈值应由漏诊与误报后果、资源、可接受工作量和公平性决定。

classification_metrics <- function(y, probability, threshold) {
  predicted <- as.integer(probability >= threshold)
  tp <- sum(predicted == 1 & y == 1)
  fp <- sum(predicted == 1 & y == 0)
  tn <- sum(predicted == 0 & y == 0)
  fn <- sum(predicted == 0 & y == 1)
  data.frame(
    阈值 = threshold,
    真阳性 = tp,
    假阳性 = fp,
    真阴性 = tn,
    假阴性 = fn,
    准确率 = (tp + tn) / (tp + fp + tn + fn),
    灵敏度 = tp / (tp + fn),
    特异度 = tn / (tn + fp),
    阳性预测值 = tp / (tp + fp)
  )
}

threshold_results <- do.call(
  rbind,
  lapply(c(0.20, 0.30, 0.50), function(threshold) {
    classification_metrics(
      test_data$infection_num,
      test_probability,
      threshold
    )
  })
)
knitr::kable(
  threshold_results,
  digits = 3,
  caption = "测试集上不同分类阈值的后果"
)
测试集上不同分类阈值的后果
阈值 真阳性 假阳性 真阴性 假阴性 准确率 灵敏度 特异度 阳性预测值
0.2 42 107 181 31 0.618 0.575 0.628 0.282
0.3 21 43 245 52 0.737 0.288 0.851 0.328
0.5 6 6 282 67 0.798 0.082 0.979 0.500

降低阈值通常提高灵敏度，同时降低特异度。准确率还可能被多数类别主导；在罕见结局中，把所有人预测为非事件也可能得到看似很高的准确率。

模型指标不是部署许可 真正应用前还需外部验证、数据漂移监测、亚组性能与校准检查、可操作性评估，以及对假阳性和假阴性后果的共同决策。受保护特征不进入模型，也不保证结果公平。

8 缺失数据、混杂与可复现性

8.1 完整案例分析会改变目标人群

glm() 默认排除模型变量中有缺失的行。只有当缺失机制和完整案例所代表的人群可以合理辩护时，这才可能合适。

set.seed(20260812)
logistic_data_missing <- train_data
missing_rows <- sample(seq_len(nrow(logistic_data_missing)), 36)
logistic_data_missing$bmi5[missing_rows] <- NA

missing_summary <- data.frame(
  总行数 = nrow(logistic_data_missing),
  BMI缺失数 = sum(is.na(logistic_data_missing$bmi5)),
  完整案例数 = sum(complete.cases(
    logistic_data_missing[c(
      "infection_num", "age10", "bmi5", "smoking", "vaccinated", "area"
    )]
  ))
)
knitr::kable(missing_summary, caption = "人为加入缺失后的样本构成")
人为加入缺失后的样本构成
总行数 BMI缺失数 完整案例数
839 36 803

应报告每个变量的缺失数量、模型实际使用人数、排除前后人群差异，并根据缺失机制考虑多重插补或敏感性分析。不要用结局均值简单填补预测变量。

8.2 变量选择不能只看 p 值

用于因果解释时，应依据因果图、时间顺序和领域知识选择混杂变量，避免调整中介或碰撞变量。用于预测时，应依据预先规定的候选变量、正则化和重采样验证控制过拟合。逐步回归会忽略模型选择带来的不确定性，并可能产生过于乐观的 p 值和性能。

8.3 一份最小可复现记录

analysis_record <- list(
  outcome_definition = "一年内呼吸道感染：1=是，0=否",
  model_formula = formula(adjusted_model),
  training_n = nobs(adjusted_model),
  event_count = sum(model.response(model.frame(adjusted_model))),
  factor_levels = lapply(
    train_data[c("smoking", "vaccinated", "area")],
    levels
  ),
  seed_for_split = 20260811,
  r_version = R.version.string
)
analysis_record
## $outcome_definition
## [1] "一年内呼吸道感染：1=是，0=否"
## 
## $model_formula
## infection_num ~ age10 + bmi5 + smoking + vaccinated + area
## <environment: 0xb0d9ef070>
## 
## $training_n
## [1] 839
## 
## $event_count
## [1] 169
## 
## $factor_levels
## $factor_levels$smoking
## [1] "否" "是"
## 
## $factor_levels$vaccinated
## [1] "否" "是"
## 
## $factor_levels$area
## [1] "城市" "农村"
## 
## 
## $seed_for_split
## [1] 20260811
## 
## $r_version
## [1] "R version 4.6.1 (2026-06-24)"

9 结果报告模板

一个清晰的报告至少回答以下问题：

  1. 人群与结局： 谁被纳入？事件如何定义？随访窗口多长？
  2. 模型： 使用什么链接函数？连续变量如何缩放或变换？分类变量参照组是什么？
  3. 调整集： 为什么选择这些变量？有多少行因缺失被排除？
  4. 估计结果： 报告 OR、95% 置信区间和明确单位；最好同时报告有意义的预测概率或绝对差。
  5. 诊断： 是否检查非线性、稀疏单元、分离、影响点和相关观测？
  6. 验证： 性能是在训练集、内部测试集还是外部数据中获得？同时报告校准与判别。
  7. 限制： 哪些偏倚、未测混杂、测量误差和可推广性问题仍存在？

四句话示例

在模拟队列中，我们使用二项分布、logit 链接的多变量逻辑回归研究一年内感染，并预先纳入年龄（每 10 年）、BMI（每 5 kg/m²）、吸烟、接种和地区。控制模型中其他变量后，接种者相对于未接种者的感染优势比为 0.52。基于同一调整模型，训练人群标准化后的预测感染概率分别为 25% 与 15.8%，但这一差值不自动具有因果含义。模型在保留测试集上的 AUC 为 0.65、Brier 分数为 0.156；仍需外部验证并评估未测混杂、数据漂移和亚组校准。

10 常见错误速查

常见说法或做法 问题 更好的做法
“OR 0.70 表示风险降低 30%” 把优势误当风险 称为优势降低 30%，并报告预测风险或风险比
直接解释 glm() 系数 默认系数在 log-odds 尺度 用 exp(coef) 得 OR，或转成概率
连续变量 OR 不写单位 读者不知道比较幅度 写明每 1、5 或 10 单位
用 0.5 作为默认最佳阈值 忽略代价与基线率 按决策后果选择并报告多阈值
只报告准确率或 AUC 忽略校准和类别不平衡 同时报概率误差、校准和阈值指标
在训练集评价性能 结果通常过于乐观 使用重采样、测试集或外部验证
逐步选择显著变量 估计与不确定性不稳定 依据问题预先指定或使用验证过的正则化
发现影响点就删除 阈值不是删除规则 回查数据，解释来源，做敏感性分析
把残差正态作为要求 套用了线性回归假设 检查 logit 函数形式、分离、影响与独立性

11 练习与答案

练习 1：从概率到优势

某人预测感染概率为 0.25。对应优势和 log-odds 是多少？

答案： 优势为 (0.25/(1-0.25)=1/3)；log-odds 为 log⁡(1/3)≈−1.10\log(1/3)\approx-1.10。
练习 2：解释连续变量 OR

模型中 age10 的 OR 为 1.40。应该怎样解释？

答案： 在模型中其他变量相同的观测之间，年龄每增加 10 年，事件优势乘以 1.40，即高 40%。这不是说事件概率增加 40%，也不自动表示年龄的因果作用。
练习 3：概率输出

为什么 predict(adjusted_model, newdata = x) 不能直接当作概率？

答案： glm 的默认预测位于链接尺度，即 log-odds。应使用 type = "response"，或对链接尺度结果应用 plogis()。
练习 4：阈值选择

若漏掉真正高风险者的代价远高于误报，阈值一般应向哪个方向移动？

答案： 通常应降低阈值以提高灵敏度，但会增加假阳性并降低特异度。最终阈值需结合资源、后续干预风险及公平性确定。
练习 5：交互项

含 vaccinated * area 时，能否把 vaccinated是 直接解释为所有地区的平均接种 OR？

答案： 不能。该系数表示参照地区（城市）中的接种 OR；农村中的比较还要加入 vaccinated是:area农村。若想要总体平均效果，应在明确定义的目标人群中计算标准化预测。
练习 6：诊断分离

某个罕见暴露组的 18 人全部未发生事件，模型给出极大的负系数和标准误。最可能的问题是什么？

答案： 可能存在完全或准完全分离。应先检查交叉表和数据质量，再考虑减少参数、增加信息或使用惩罚/偏倚校正方法；不能把巨大 OR 当成稳定证据。

12 快速参考

目标 基础 R 写法
拟合模型 glm(y ~ x1 + x2, family = binomial(), data = d)
查看 log-odds 系数 coef(model)
计算 OR exp(coef(model))
Wald OR 置信区间 exp(coef(model) + outer(sqrt(diag(vcov(model))), qnorm(c(.025, .975))))（注意整理维度）
预测概率 predict(model, newdata = d_new, type = "response")
链接尺度与标准误 predict(model, newdata = d_new, type = "link", se.fit = TRUE)
偏差残差 residuals(model, type = "deviance")
杠杆值 hatvalues(model)
Cook 距离 cooks.distance(model)
嵌套模型比较 anova(model_small, model_large, test = "Chisq")
信息准则 AIC(model)

最终检查清单

本教程提供逻辑回归的学习框架，而不是临床或公共卫生决策规则。真实分析应结合研究设计、抽样结构、领域知识、预注册计划和适当的软件验证。