The V Lab
适用对象医学、护理、公共卫生与临床研究学习者
学习时长约 150–210 分钟
先修要求描述统计、置信区间与基础 R

关于数据与用途 本教程中的所有个体记录、效应和检验结果均由固定随机种子模拟，只用于统计学教学，不包含真实患者信息，也不能作为任何治疗有效或安全的临床证据。

如何使用本教程

本教程不是一张按变量名称机械查找检验的菜单。每一节都遵循同一条路径：明确研究问题与目标估计量 → 识别研究设计和数据结构 → 检查数据质量与方法假设 → 估计效应及其不确定性 → 必要时进行检验 → 结合临床意义报告。

建议先阅读“检验选择总览”，再学习与你的数据结构相关的章节。代码默认显示，可逐块运行；每个例子都尽量同时给出效应量、95% 置信区间和 p 值。

学习目标

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

  • 根据结局类型、组数、独立或配对结构以及研究设计选择常见统计方法；
  • 区分单样本、两独立样本、配对样本、多组独立样本和重复测量问题；
  • 正确使用并解释 t 检验、Wilcoxon 检验、方差分析和 Kruskal–Wallis 检验；
  • 正确使用卡方检验、Fisher 精确检验、McNemar 检验和趋势检验；
  • 根据问题选择 Pearson、Spearman 或 Kendall 相关；
  • 使用 log-rank 检验比较生存曲线，并理解其不能替代效应估计；
  • 检查关键假设，而不把“先做正态性检验”当作自动决策器；
  • 处理多重比较，区分统计显著、临床重要、等效和非劣效；
  • 报告效应大小、方向、95% 置信区间、p 值、样本量、假设与局限。

1 先从研究问题选择方法

1.1 统计检验不是研究设计的补丁

统计检验衡量数据与某个零假设及其模型假设的相容程度。它不能修复选择偏倚、失访、错误测量、未控制混杂、错误的时间顺序或不合适的研究问题。

在计算任何 p 值前，应写清楚：

  1. 总体与分析单位：每行代表患者、眼睛、住院、病灶还是一次访视？
  2. 结局：连续值、二分类、多分类、有序等级、计数还是事件时间？
  3. 比较结构：一组与参考值、两个独立组、同一患者前后、多组还是重复测量？
  4. 目标估计量：均值差、风险差、风险比、优势比、相关系数还是生存概率差？
  5. 依赖结构：观测是否配对、聚类、重复、删失或经过抽样加权？
  6. 临床阈值：多大的差异才具有临床或公共卫生意义？

最常见的起点错误 “我的变量正态吗，所以该用哪种检验？”通常不是第一个问题。应先确认比较对象、独立性、配对关系、结局尺度和目标估计量；这些因素往往比边际分布是否完美正态更重要。

1.2 快速选择矩阵

研究问题与结构 常用主要方法 主要效应或估计量 重要提醒
一组连续值与参考值比较 单样本 t 检验 均值与参考值之差 对均值推断；检查异常值与抽样机制
两个独立组的连续结局 Welch t 检验 均值差 默认不要求两组方差相等
同一对象前后连续结局 配对 t 检验 个体差值的均值 分析的是“差值”，不能当成独立组
两独立组的有序或明显偏态结局 Wilcoxon 秩和检验 分布位置／概率优势 通常不是“中位数检验”
配对有序或偏态差值 Wilcoxon 符号秩检验 差值分布的位置 需要差值分布近似对称；零值处理要说明
三个及以上独立组连续结局 单因素／Welch ANOVA 组均值差异 总体检验后再做预设或校正的比较
三个及以上独立组有序或偏态结局 Kruskal–Wallis 检验 秩分布差异 显著后需校正的两两比较
三个及以上重复测量 重复测量模型／Friedman 检验 时间或条件差异 必须保留患者内相关性
两个独立分类变量 Pearson 卡方检验 比例关联 小期望频数时考虑 Fisher 精确检验
两个独立组的稀疏 2×2 表 Fisher 精确检验 优势比及精确推断 不等于风险比；设计固定边际的假设要理解
同一对象前后二分类结局 McNemar 检验 不一致配对的方向差 只使用不一致配对的信息
两个连续变量的线性关联 Pearson 相关 相关系数 rr 检查散点图、异常值与独立性
单调但非线性／有序关联 Spearman 或 Kendall 相关 秩相关 不是因果效应，也不等于一致性
两组删失事件时间 Kaplan–Meier + log-rank 生存曲线差异 同时报告生存概率、时间和风险集

核心原则：检验回答“数据与零假设是否相容”，效应估计回答“差异有多大、方向如何”。医学报告通常需要后者及其置信区间，而不仅是一颗显著性星号。

1.3 独立、配对与聚类必须先分清

  • 独立样本：两组由不同患者组成，一位患者只贡献一个独立观察。
  • 配对样本：同一患者前后、左右眼、匹配病例与对照等观测具有明确对应关系。
  • 重复测量：一位患者在三个或更多时点贡献数据。
  • 聚类数据：患者嵌套于医院、牙齿嵌套于患者、多个病灶嵌套于同一人。

普通 t 检验、卡方检验和秩检验通常假定分析单位相互独立。把同一患者的多个观测当成独立样本会低估标准误并夸大证据。复杂重复或聚类资料一般需要混合效应模型、广义估计方程、聚类稳健标准误或设计相匹配的随机化方法。

2 连续结局：一组与两组比较

2.1 单样本 t 检验

单样本 t 检验用于比较一个总体的均值与预先指定的参考值 μ0\mu_0：

H0:μ=μ0,t=x‾−μ0s/n. H_0:\mu=\mu_0,\qquad t=\frac{\bar x-\mu_0}{s/\sqrt{n}}.

示例问题：这组患者的平均 12 周收缩压下降是否不同于预先设定的 8 mmHg？参考值应来自方案、历史标准或临床问题，不能在看完数据后挑选。

new_treatment_reduction <- subset(
  trial_data,
  treatment == "新治疗"
)$sbp_reduction

one_sample_result <- t.test(
  new_treatment_reduction,
  mu = 8,
  conf.level = 0.95
)

data.frame(
  样本量 = length(new_treatment_reduction),
  样本均值 = mean(new_treatment_reduction),
  均值与8的差 = mean(new_treatment_reduction) - 8,
  CI下限 = unname(one_sample_result$conf.int[1] - 8),
  CI上限 = unname(one_sample_result$conf.int[2] - 8),
  p值 = one_sample_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(
    digits = 3,
    caption = "新治疗组平均收缩压下降与 8 mmHg 的比较"
  )
新治疗组平均收缩压下降与 8 mmHg 的比较
样本量 样本均值 均值与8的差 CI下限 CI上限 p值
160 11.6 3.62 2.23 5.01 0

关键条件包括：观察独立；结局的均值有明确意义；没有足以支配均值和标准误的严重错误值；样本来自与推断目标相符的过程。小样本时还需差值分布近似正态；样本量较大时，均值的抽样分布通常更稳健，但极端偏态和异常值仍需认真处理。

2.2 两独立样本：优先考虑 Welch t 检验

当研究目标是两个独立组的均值差时，Welch t 检验通常是合理默认选择。它不强迫两组总体方差相等，在方差确实相等时效率损失通常很小。

welch_result <- t.test(
  sbp_reduction ~ treatment,
  data = trial_data,
  var.equal = FALSE
)

group_summary <- do.call(
  rbind,
  lapply(split(trial_data$sbp_reduction, trial_data$treatment), function(x) {
    c(n = length(x), mean = mean(x), sd = sd(x))
  })
)
group_summary <- data.frame(
  treatment = rownames(group_summary),
  group_summary,
  row.names = NULL
)

knitr::kable(
  group_summary,
  digits = 2,
  col.names = c("治疗组", "n", "均值", "标准差"),
  caption = "各组收缩压下降的描述统计"
)
各组收缩压下降的描述统计
治疗组 n 均值 标准差
标准治疗 160 7.07 8.24
新治疗 160 11.62 8.90 
data.frame(
  对比 = "标准治疗 − 新治疗",
  均值差 = unname(welch_result$estimate[1] - welch_result$estimate[2]),
  CI下限 = unname(welch_result$conf.int[1]),
  CI上限 = unname(welch_result$conf.int[2]),
  自由度 = unname(welch_result$parameter),
  p值 = welch_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(
    digits = 3,
    caption = "Welch 两独立样本 t 检验"
  )
Welch 两独立样本 t 检验
对比 均值差 CI下限 CI上限 自由度 p值
标准治疗 − 新治疗 -4.55 -6.43 -2.66 316 0

t.test() 中因子第一水平减去第二水平，所以这里的差值是“标准治疗 − 新治疗”。如果希望正值代表新治疗获得更大下降，应在报告表中明确反转方向，不能只抄输出而不核对编码。

2.2.1 效应量：均值差优先，标准化差异辅助

原始单位的均值差通常最容易判断临床意义。标准化均值差可用于不同量表间比较，但会受到研究人群变异程度影响。

standard_values <- subset(
  trial_data,
  treatment == "标准治疗"
)$sbp_reduction
new_values <- subset(
  trial_data,
  treatment == "新治疗"
)$sbp_reduction

n0 <- length(standard_values)
n1 <- length(new_values)
pooled_sd <- sqrt(
  ((n0 - 1) * var(standard_values) + (n1 - 1) * var(new_values)) /
    (n0 + n1 - 2)
)
cohens_d <- (mean(new_values) - mean(standard_values)) / pooled_sd
small_sample_correction <- 1 - 3 / (4 * (n0 + n1) - 9)
hedges_g <- small_sample_correction * cohens_d

data.frame(
  指标 = c("新治疗 − 标准治疗的均值差（mmHg）", "Hedges g"),
  估计值 = c(mean(new_values) - mean(standard_values), hedges_g)
) |>
  knitr::kable(digits = 3, caption = "原始与标准化效应量")
原始与标准化效应量
指标 估计值
新治疗 − 标准治疗的均值差（mmHg） 4.548
Hedges g 0.529

推荐报告句式 新治疗组与标准治疗组各纳入多少人；分别报告均值（SD）；随后报告预先规定方向的均值差、95% CI 与 Welch t 检验 p 值；最后说明该区间相对于最小临床重要差异意味着什么。

2.3 配对 t 检验

配对 t 检验将每位患者的“治疗后 − 治疗前”先变为一个差值，再检验平均差值是否为零。它假定不同患者的差值相互独立，并关注差值分布，而不是分别要求治疗前和治疗后都正态。

paired_result <- t.test(
  trial_data$followup_sbp,
  trial_data$baseline_sbp,
  paired = TRUE
)

within_person_change <- trial_data$followup_sbp - trial_data$baseline_sbp

data.frame(
  样本量 = length(within_person_change),
  平均变化 = mean(within_person_change),
  变化SD = sd(within_person_change),
  CI下限 = unname(paired_result$conf.int[1]),
  CI上限 = unname(paired_result$conf.int[2]),
  p值 = paired_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(
    digits = 3,
    caption = "全体参与者随访值减基线值的配对 t 检验"
  )
全体参与者随访值减基线值的配对 t 检验
样本量 平均变化 变化SD CI下限 CI上限 p值
320 -9.35 8.86 -10.3 -8.37 0 
plot_ids <- seq_len(50)
matplot(
  x = c(1, 2),
  y = t(as.matrix(trial_data[plot_ids, c("baseline_sbp", "followup_sbp")])),
  type = "l",
  lty = 1,
  col = grDevices::adjustcolor(palette_test[["gray"]], alpha.f = 0.32),
  xaxt = "n",
  xlab = "时点",
  ylab = "收缩压（mmHg）"
)
axis(1, at = c(1, 2), labels = c("基线", "第 12 周"))
points(
  c(1, 2),
  c(mean(trial_data$baseline_sbp), mean(trial_data$followup_sbp)),
  type = "b",
  pch = 19,
  lwd = 3,
  col = palette_test[["vermillion"]]
)
配对折线图，显示每位参与者从基线到十二周的收缩压变化。

每位参与者的基线与 12 周收缩压。细线表示同一人。

前后显著性不同不等于组间变化不同 在一个治疗组内“前后显著”、另一个组内“前后不显著”，不能证明两组变化不同。应直接比较两组的个体变化，或更常见地用基线调整模型比较随访结局，并保留随机化分组。

2.4 Wilcoxon 秩和检验：两个独立组

当结局是有序变量，或连续结局存在强烈偏态而秩比较与问题相符时，可使用 Wilcoxon 秩和检验（Mann–Whitney U 检验）。它利用两组混合排序后的秩。

rank_sum_result <- wilcox.test(
  crp_week12 ~ treatment,
  data = trial_data,
  exact = FALSE,
  conf.int = TRUE
)

crp_summary <- do.call(
  rbind,
  lapply(split(trial_data$crp_week12, trial_data$treatment), function(x) {
    c(
      n = length(x),
      median = median(x),
      q1 = unname(quantile(x, 0.25)),
      q3 = unname(quantile(x, 0.75))
    )
  })
)

knitr::kable(
  data.frame(治疗组 = rownames(crp_summary), crp_summary, row.names = NULL),
  digits = 2,
  col.names = c("治疗组", "n", "中位数", "Q1", "Q3"),
  caption = "12 周 CRP 的中位数与四分位数"
)
12 周 CRP 的中位数与四分位数
治疗组 n 中位数 Q1 Q3
标准治疗 160 3.09 1.96 5.53
新治疗 160 2.88 1.62 4.69 
data.frame(
  W统计量 = unname(rank_sum_result$statistic),
  位置差估计 = unname(rank_sum_result$estimate),
  CI下限 = unname(rank_sum_result$conf.int[1]),
  CI上限 = unname(rank_sum_result$conf.int[2]),
  p值 = rank_sum_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(
    digits = 3,
    caption = "Wilcoxon 秩和检验及 Hodges–Lehmann 型位置差"
  )
Wilcoxon 秩和检验及 Hodges–Lehmann 型位置差
W统计量 位置差估计 CI下限 CI上限 p值
14232 0.36 -0.06 0.81 0.084

若两组分布形状相同、主要差别可由整体平移描述，结果可解释为位置差。一般情况下，它检验的是两组分布是否具有相同的秩结构，不能自动写成中位数相等的检验。即使使用秩检验，也要画图并报告各组分布。

boxplot(
  crp_week12 ~ treatment,
  data = trial_data,
  col = c(palette_test[["sky"]], palette_test[["orange"]]),
  border = palette_test[["navy"]],
  ylab = "CRP（mg/L）",
  xlab = "治疗组",
  outline = FALSE
)
stripchart(
  crp_week12 ~ treatment,
  data = trial_data,
  vertical = TRUE,
  method = "jitter",
  pch = 16,
  cex = 0.55,
  col = grDevices::adjustcolor(palette_test[["navy"]], alpha.f = 0.32),
  add = TRUE
)
按治疗组分面的箱线图和抖动散点，显示十二周 CRP 的右偏分布。

两组 12 周 CRP 的右偏分布。

2.5 Wilcoxon 符号秩与符号检验：配对资料

Wilcoxon 符号秩检验适用于配对差值的秩分析。它除独立配对外，还通常需要差值分布围绕中心近似对称。若对称性难以辩护，可用只看正负号的精确符号检验，但它会丢弃差值大小信息。

signed_rank_result <- wilcox.test(
  trial_data$followup_sbp,
  trial_data$baseline_sbp,
  paired = TRUE,
  exact = FALSE,
  conf.int = TRUE
)

nonzero_change <- within_person_change[within_person_change != 0]
negative_changes <- sum(nonzero_change < 0)
sign_test_result <- binom.test(
  negative_changes,
  length(nonzero_change),
  p = 0.5
)

data.frame(
  方法 = c("Wilcoxon 符号秩检验", "精确符号检验"),
  统计量或计数 = c(
    unname(signed_rank_result$statistic),
    negative_changes
  ),
  p值 = c(signed_rank_result$p.value, sign_test_result$p.value),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "两种配对非参数检验")
两种配对非参数检验
方法 统计量或计数 p值
Wilcoxon 符号秩检验 3260 0
精确符号检验 276 0
检查理解：t 检验还是秩检验？

一项随机试验比较两组患者住院天数。大量患者住院 1–3 天，少数患者住院超过 60 天。研究者最关心平均住院资源占用。应否因为分布偏态就自动改用 Wilcoxon 检验？

答案：不应自动更换。 若目标估计量是平均住院天数，Wilcoxon 检验回答的不是同一个问题。可先核对极端值是否真实，再考虑稳健或自助法均值差区间、置换方法、适合正偏结局的模型，以及同时报告均值与分布。方法必须与目标估计量一致。

3 三组及以上的连续或有序结局

3.1 单因素 ANOVA 与 Welch ANOVA

单因素方差分析的总体零假设是所有组的总体均值相等。它回答“是否至少有一个组均值不同”，但不会自动告诉我们具体哪些组不同。传统 ANOVA 假定独立误差、合理的均值模型、组内残差近似正态和总体方差相等；Welch ANOVA 放宽方差相等要求。

dose_summary <- do.call(
  rbind,
  lapply(split(dose_data$hemoglobin_change, dose_data$dose), function(x) {
    c(n = length(x), mean = mean(x), sd = sd(x))
  })
)

knitr::kable(
  data.frame(组别 = rownames(dose_summary), dose_summary, row.names = NULL),
  digits = 3,
  col.names = c("组别", "n", "均值", "标准差"),
  caption = "三组血红蛋白变化的描述统计"
)
三组血红蛋白变化的描述统计
组别 n 均值 标准差
安慰剂 60 0.075 0.773
低剂量 60 0.695 0.847
高剂量 60 1.219 1.146 
traditional_anova <- aov(hemoglobin_change ~ dose, data = dose_data)
welch_anova <- oneway.test(
  hemoglobin_change ~ dose,
  data = dose_data,
  var.equal = FALSE
)

anova_table <- summary(traditional_anova)[[1]]
ss_between <- anova_table["dose", "Sum Sq"]
ss_total <- sum(anova_table[, "Sum Sq"])
ms_error <- anova_table["Residuals", "Mean Sq"]
k_groups <- nlevels(dose_data$dose)
omega_squared <- max(
  0,
  (ss_between - (k_groups - 1) * ms_error) / (ss_total + ms_error)
)

data.frame(
  方法 = c("传统单因素 ANOVA", "Welch ANOVA"),
  统计量 = c(
    anova_table["dose", "F value"],
    unname(welch_anova$statistic)
  ),
  分子自由度 = c(
    anova_table["dose", "Df"],
    unname(welch_anova$parameter[1])
  ),
  分母自由度 = c(
    anova_table["Residuals", "Df"],
    unname(welch_anova$parameter[2])
  ),
  p值 = c(
    anova_table["dose", "Pr(>F)"],
    welch_anova$p.value
  ),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "传统 ANOVA 与 Welch ANOVA")
传统 ANOVA 与 Welch ANOVA
方法 统计量 分子自由度 分母自由度 p值
传统单因素 ANOVA 22.5 2 177 0
Welch ANOVA 22.3 2 115 0 
data.frame(
  效应量 = "omega-squared",
  估计值 = omega_squared
) |>
  knitr::kable(digits = 3, caption = "传统单因素 ANOVA 的总体效应量")
传统单因素 ANOVA 的总体效应量
效应量 估计值
omega-squared 0.193

ω2\omega^2 近似描述结局总变异中可归因于组别差异的比例，但不应脱离研究背景套用固定“小、中、大”阈值。它也不能代替具体组间均值差。

3.1.1 总体检验后的比较

若方案预先指定“低剂量对安慰剂”和“高剂量对安慰剂”，这些对比比先做一次总体检验、显著后再穷举所有组合更贴近研究问题。若确实关心全部两两比较，可使用 Tukey 方法；方差不等时可进行各组 Welch t 检验并校正多重性。

tukey_result <- TukeyHSD(traditional_anova, "dose")$dose
knitr::kable(
  data.frame(比较 = rownames(tukey_result), tukey_result, row.names = NULL),
  digits = 3,
  col.names = c("比较", "均值差", "同时CI下限", "同时CI上限", "调整后p值"),
  caption = "Tukey 全部两两比较"
)
Tukey 全部两两比较
比较 均值差 同时CI下限 同时CI上限 调整后p值
低剂量-安慰剂 0.620 0.216 1.024 0.001
高剂量-安慰剂 1.144 0.740 1.548 0.000
高剂量-低剂量 0.524 0.120 0.928 0.007 
pairwise_welch <- pairwise.t.test(
  dose_data$hemoglobin_change,
  dose_data$dose,
  p.adjust.method = "holm",
  pool.sd = FALSE
)
pairwise_welch$p.value |>
  knitr::kable(
    digits = 3,
    caption = "两两 Welch t 检验的 Holm 调整 p 值"
  )
两两 Welch t 检验的 Holm 调整 p 值
安慰剂 低剂量
低剂量 0 NA
高剂量 0 0.005

注意：pairwise.t.test() 上表提供调整后的 p 值，但没有同时置信区间。正式报告若以多重对比为主要结论，应使用能生成与校正方法一致的同时区间的模型后对比工具，而不是并列展示调整 p 值与未调整区间。

3.2 Kruskal–Wallis 检验

Kruskal–Wallis 检验将 Wilcoxon 秩和方法扩展到三个及以上独立组。其零假设涉及各组秩分布相同；只有在分布形状相近、差异可视为位置平移时，才适合简化为位置比较。

kruskal_result <- kruskal.test(
  inflammation_score ~ dose,
  data = dose_data
)

h_statistic <- unname(kruskal_result$statistic)
epsilon_squared <- max(
  0,
  (h_statistic - k_groups + 1) / (nrow(dose_data) - k_groups)
)

rank_summary <- do.call(
  rbind,
  lapply(split(dose_data$inflammation_score, dose_data$dose), function(x) {
    c(
      n = length(x),
      median = median(x),
      q1 = unname(quantile(x, 0.25)),
      q3 = unname(quantile(x, 0.75))
    )
  })
)

knitr::kable(
  data.frame(组别 = rownames(rank_summary), rank_summary, row.names = NULL),
  digits = 2,
  col.names = c("组别", "n", "中位数", "Q1", "Q3"),
  caption = "三组炎症评分的描述统计"
)
三组炎症评分的描述统计
组别 n 中位数 Q1 Q3
安慰剂 60 2.96 1.40 5.43
低剂量 60 2.51 1.65 3.70
高剂量 60 2.58 1.60 4.14 
data.frame(
  H统计量 = h_statistic,
  自由度 = unname(kruskal_result$parameter),
  p值 = kruskal_result$p.value,
  秩epsilon平方 = epsilon_squared,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Kruskal–Wallis 检验")
Kruskal–Wallis 检验
H统计量 自由度 p值 秩epsilon平方
1.46 2 0.483 0 
pairwise_rank <- pairwise.wilcox.test(
  dose_data$inflammation_score,
  dose_data$dose,
  p.adjust.method = "holm",
  exact = FALSE
)

pairwise_rank$p.value |>
  knitr::kable(
    digits = 3,
    caption = "两两 Wilcoxon 秩和检验的 Holm 调整 p 值"
  )
两两 Wilcoxon 秩和检验的 Holm 调整 p 值
安慰剂 低剂量
低剂量 0.839 NA
高剂量 0.839 0.967

总体显著不是每一对都不同 ANOVA 或 Kruskal–Wallis 的小 p 值只说明数据不支持“所有组完全相同”的总体零假设。具体组间差异必须通过预先指定或经过恰当多重性控制的对比来估计。

3.3 三个及以上重复测量

同一患者在多个时点的观测通常正相关。完整且平衡的连续资料可用重复测量 ANOVA；Friedman 检验是完整区组的秩方法。真实纵向研究若存在不等时间、缺失随访、时变协变量或组间比较，混合效应模型或广义估计方程通常更合适。

repeated_summary <- do.call(
  rbind,
  lapply(split(repeated_data$pain_score, repeated_data$time), function(x) {
    c(n = length(x), mean = mean(x), sd = sd(x), median = median(x))
  })
)

knitr::kable(
  data.frame(时点 = rownames(repeated_summary), repeated_summary, row.names = NULL),
  digits = 2,
  col.names = c("时点", "n", "均值", "标准差", "中位数"),
  caption = "三个时点疼痛评分的描述统计"
)
三个时点疼痛评分的描述统计
时点 n 均值 标准差 中位数
基线 90 6.53 1.59 6.55
第4周 90 4.97 1.59 5.05
第12周 90 4.08 1.48 4.10 
repeated_anova <- aov(
  pain_score ~ time + Error(patient_id / time),
  data = repeated_data
)
summary(repeated_anova)
## 
## Error: patient_id
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 89    541    6.08               
## 
## Error: patient_id:time
##            Df Sum Sq Mean Sq F value Pr(>F)    
## time        2    276   137.9     235 <2e-16 ***
## Residuals 178    104     0.6                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

上述经典重复测量 ANOVA 对三个以上时点还涉及球形性假设。aov() 的简化输出不自动提供 Mauchly 检验或 Greenhouse–Geisser 校正，因此在需要正式推断时应使用能明确处理球形性或协方差结构的方法。

friedman_result <- friedman.test(
  pain_score ~ time | patient_id,
  data = repeated_data
)

data.frame(
  卡方统计量 = unname(friedman_result$statistic),
  自由度 = unname(friedman_result$parameter),
  p值 = friedman_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Friedman 重复测量秩检验")
Friedman 重复测量秩检验
卡方统计量 自由度 p值
134 2 0

如需探索哪些时点不同，可对配对比较的 p 值进行校正：

time_pairs <- combn(colnames(pain_matrix), 2, simplify = FALSE)
raw_p <- vapply(time_pairs, function(pair) {
  wilcox.test(
    pain_matrix[, pair[1]],
    pain_matrix[, pair[2]],
    paired = TRUE,
    exact = FALSE
  )$p.value
}, numeric(1))

repeated_pairwise <- data.frame(
  比较 = vapply(time_pairs, paste, collapse = " vs ", FUN.VALUE = character(1)),
  原始p值 = raw_p,
  Holm调整p值 = p.adjust(raw_p, method = "holm"),
  check.names = FALSE
)

knitr::kable(
  repeated_pairwise,
  digits = 3,
  caption = "三个时点配对秩检验的多重性校正"
)
三个时点配对秩检验的多重性校正
比较 原始p值 Holm调整p值
基线 vs 第4周 0 0
基线 vs 第12周 0 0
第4周 vs 第12周 0 0

随机试验有基线测量时 两组随机试验通常不应只比较各组内部前后 p 值。对于连续随访结局，预先规定的基线调整 ANCOVA 常能直接比较组间随访均值、提高精度并处理偶然基线不平衡；若有多个随访时点，再考虑纵向模型。

4 分类结局与比例

4.1 单个比例：精确二项与得分近似

若 80 名患者中 62 名达到治疗反应，可比较真实反应比例是否为预设值 0.70。binom.test() 使用精确二项分布，prop.test() 给出大样本得分型近似；两者的区间和 p 值可能略有不同。

x_response <- 62
n_response <- 80
exact_proportion <- binom.test(x_response, n_response, p = 0.70)
score_proportion <- prop.test(
  x_response,
  n_response,
  p = 0.70,
  correct = FALSE
)

data.frame(
  方法 = c("精确二项", "得分近似"),
  观察比例 = x_response / n_response,
  CI下限 = c(exact_proportion$conf.int[1], score_proportion$conf.int[1]),
  CI上限 = c(exact_proportion$conf.int[2], score_proportion$conf.int[2]),
  p值 = c(exact_proportion$p.value, score_proportion$p.value),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "单个反应比例的两种推断方法")
单个反应比例的两种推断方法
方法 观察比例 CI下限 CI上限 p值
精确二项 0.775 0.668 0.861 0.179
得分近似 0.775 0.672 0.853 0.143

4.2 两个独立比例：卡方、比例检验与效应量

二分类结局不能只报告卡方 p 值。随机试验或队列研究中，风险差（RD）和风险比（RR）通常更直接；病例对照研究由于按结局抽样，通常以优势比（OR）为主要关联尺度。

response_table <- table(trial_data$treatment, trial_data$response)
response_table |>
  addmargins() |>
  knitr::kable(caption = "治疗组与 12 周反应的 2×2 表")
治疗组与 12 周反应的 2×2 表
否 是 Sum
标准治疗 67 93 160
新治疗 49 111 160
Sum 116 204 320 
response_risks <- prop.table(response_table, margin = 1)[, "是"]
data.frame(
  治疗组 = names(response_risks),
  反应风险 = unname(response_risks)
) |>
  knitr::kable(digits = 3, caption = "各组观察到的反应风险")
各组观察到的反应风险
治疗组 反应风险
标准治疗 0.581
新治疗 0.694 
chi_result <- chisq.test(response_table, correct = FALSE)
prop_result <- prop.test(
  x = response_table[, "是"],
  n = rowSums(response_table),
  correct = FALSE
)

data.frame(
  方法 = c("Pearson 卡方检验", "两样本比例检验"),
  统计量 = c(unname(chi_result$statistic), unname(prop_result$statistic)),
  自由度 = c(unname(chi_result$parameter), unname(prop_result$parameter)),
  p值 = c(chi_result$p.value, prop_result$p.value),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "两独立比例的总体检验")
两独立比例的总体检验
方法 统计量 自由度 p值
Pearson 卡方检验 4.38 1 0.036
两样本比例检验 4.38 1 0.036 
risk_measures(response_table) |>
  knitr::kable(
    digits = 3,
    caption = "反应结局的绝对与相对效应（大样本 Wald 区间）"
  )
反应结局的绝对与相对效应（大样本 Wald 区间）
指标 估计值 CI下限 CI上限
新治疗风险 0.694 NA NA
标准治疗风险 0.581 NA NA
风险差 0.112 0.008 0.217
风险比 1.194 1.010 1.411
优势比 1.632 1.030 2.585

风险差的零值是 0，风险比和优势比的零效应值是 1。上述手工区间用于演示；稀疏数据、聚类资料或调整后分析应使用与设计和模型匹配的方法。

4.2.1 期望频数与 Cramér’s V

卡方近似依赖期望频数足够支持渐近分布，不能用“总样本量小于 30”这样的固定口号判断。先查看期望频数和稀疏程度：

knitr::kable(
  chi_result$expected,
  digits = 2,
  caption = "独立性零假设下的期望频数"
)
独立性零假设下的期望频数
否 是
标准治疗 58 102
新治疗 58 102 
cramers_v <- sqrt(
  unname(chi_result$statistic) /
    (sum(response_table) * min(nrow(response_table) - 1, ncol(response_table) - 1))
)
data.frame(效应量 = "Cramér's V", 估计值 = cramers_v) |>
  knitr::kable(digits = 3, caption = "分类关联的标准化效应量")
分类关联的标准化效应量
效应量 估计值
Cramér’s V 0.117

4.3 Fisher 精确检验：稀疏 2×2 表

罕见严重不良事件常产生小单元格。Fisher 精确检验可避免依赖大样本卡方近似，但仍需报告原始分母、绝对风险与效应区间。

rare_event_table <- matrix(
  c(1, 39, 7, 33),
  nrow = 2,
  byrow = TRUE,
  dimnames = list(
    治疗 = c("新治疗", "标准治疗"),
    严重不良事件 = c("是", "否")
  )
)

fisher_result <- fisher.test(rare_event_table)

knitr::kable(
  addmargins(rare_event_table),
  caption = "稀疏严重不良事件 2×2 表"
)
稀疏严重不良事件 2×2 表
是 否 Sum
新治疗 1 39 40
标准治疗 7 33 40
Sum 8 72 80 
data.frame(
  条件优势比估计 = unname(fisher_result$estimate),
  CI下限 = unname(fisher_result$conf.int[1]),
  CI上限 = unname(fisher_result$conf.int[2]),
  p值 = fisher_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Fisher 精确检验")
Fisher 精确检验
条件优势比估计 CI下限 CI上限 p值
0.124 0.003 1.04 0.057

fisher.test() 返回的条件最大似然 OR 估计可能与简单交叉乘积 ad/bcad/bc 略有不同。零单元格还可能使简单 OR 为 0 或无穷；不应只为得到有限数字而随意加 0.5，除非分析方法已预先说明并有合理依据。

4.4 配对二分类结局：McNemar 检验

同一患者干预前后是否有症状是配对二分类资料。McNemar 检验只使用两类不一致配对：前有后无与前无后有。

knitr::kable(
  addmargins(paired_symptom),
  caption = "同一患者干预前后的症状状态"
)
同一患者干预前后的症状状态
无症状 有症状 Sum
无症状 49 6 55
有症状 45 40 85
Sum 94 46 140 
mcnemar_result <- mcnemar.test(paired_symptom, correct = TRUE)
b <- paired_symptom["无症状", "有症状"]
c <- paired_symptom["有症状", "无症状"]
exact_mcnemar <- binom.test(b, b + c, p = 0.5)

data.frame(
  方法 = c("McNemar 卡方近似（连续性校正）", "不一致配对的精确二项检验"),
  统计量或计数 = c(unname(mcnemar_result$statistic), b),
  p值 = c(mcnemar_result$p.value, exact_mcnemar$p.value),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "配对二分类结局的检验")
配对二分类结局的检验
方法 统计量或计数 p值
McNemar 卡方近似（连续性校正） 28.3 0
不一致配对的精确二项检验 6.0 0 
before_risk <- mean(symptom_before)
after_risk <- mean(symptom_after)
paired_risk_difference <- after_risk - before_risk
matched_or <- b / c

data.frame(
  指标 = c("干预前症状比例", "干预后症状比例", "后 − 前风险差", "不一致配对优势比"),
  估计值 = c(before_risk, after_risk, paired_risk_difference, matched_or)
) |>
  knitr::kable(digits = 3, caption = "配对二分类资料的描述性效应")
配对二分类资料的描述性效应
指标 估计值
干预前症状比例 0.607
干预后症状比例 0.329
后 − 前风险差 -0.279
不一致配对优势比 0.133

若不一致对很少，应优先报告精确结果，并诚实呈现宽区间。配对 OR 只描述不一致对方向，不能与独立 2×2 表 OR 混为一谈。

4.5 有序暴露的比例趋势

当剂量水平有科学上预先规定的顺序，prop.trend.test() 可检验比例是否呈线性趋势。它比无序卡方检验更聚焦，但不能发现所有非线性模式。

responders_by_dose <- c(21, 31, 43)
total_by_dose <- c(60, 60, 60)
trend_result <- prop.trend.test(
  responders_by_dose,
  total_by_dose,
  score = c(0, 1, 2)
)

data.frame(
  剂量 = c("安慰剂", "低剂量", "高剂量"),
  反应人数 = responders_by_dose,
  总人数 = total_by_dose,
  反应比例 = responders_by_dose / total_by_dose
) |>
  knitr::kable(digits = 3, caption = "按有序剂量分组的反应比例")
按有序剂量分组的反应比例
剂量 反应人数 总人数 反应比例
安慰剂 21 60 0.350
低剂量 31 60 0.517
高剂量 43 60 0.717 
data.frame(
  趋势卡方 = unname(trend_result$statistic),
  自由度 = unname(trend_result$parameter),
  p值 = trend_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "Cochran–Armitage 型比例趋势检验")
Cochran–Armitage 型比例趋势检验
趋势卡方 自由度 p值
16.2 1 0

4.6 分层 2×2 表：Mantel–Haenszel 检验

若治疗与结局的 2×2 表按中心或一个预先规定的混杂因素分层，可用 Mantel–Haenszel 方法估计共同 OR 并检验条件关联。它假定各层 OR 足够相似；若效应在各层明显不同，一个共同 OR 会掩盖效应修饰。

mh_array <- array(
  c(
    28, 22, 19, 31,
    34, 16, 25, 25,
    19, 31, 13, 37
  ),
  dim = c(2, 2, 3),
  dimnames = list(
    治疗 = c("新治疗", "标准治疗"),
    反应 = c("是", "否"),
    中心 = c("中心A", "中心B", "中心C")
  )
)

mh_result <- mantelhaen.test(mh_array)

data.frame(
  共同优势比 = unname(mh_result$estimate),
  CI下限 = unname(mh_result$conf.int[1]),
  CI上限 = unname(mh_result$conf.int[2]),
  p值 = mh_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "按研究中心分层的 Mantel–Haenszel 分析")
按研究中心分层的 Mantel–Haenszel 分析
共同优势比 CI下限 CI上限 p值
1.98 1.24 3.18 0.007
检查理解：卡方还是 McNemar？

一项研究让每名患者先后接受两种快速检测，并比较两种检测的阳性率。应使用 Pearson 卡方检验还是 McNemar 检验？

答案：McNemar 检验。 两种检测来自同一名患者，结果配对。普通 Pearson 卡方检验会错误地把两组结果视为独立。若研究目标是比较灵敏度，则分析还必须限于按参考标准确诊者，并继续保留配对结构。

5 事件计数与发病率

5.1 Poisson 率检验

当各参与者随访时间不同，单纯比较“发生过事件的比例”可能丢失信息。若关注事件发生速率，应报告事件数、总风险人时和明确单位下的发病率。

假设新治疗组在 410 人年中发生 18 次事件，标准治疗组在 395 人年中发生 31 次事件：

event_counts <- c(新治疗 = 18, 标准治疗 = 31)
person_years <- c(新治疗 = 410, 标准治疗 = 395)
rate_result <- poisson.test(event_counts, person_years)

rate_table <- data.frame(
  组别 = names(event_counts),
  事件数 = unname(event_counts),
  人年 = unname(person_years),
  每100人年发生率 = 100 * unname(event_counts / person_years)
)

knitr::kable(
  rate_table,
  digits = 2,
  caption = "两组事件数、人时和发生率"
)
两组事件数、人时和发生率
组别 事件数 人年 每100人年发生率
新治疗 18 410 4.39
标准治疗 31 395 7.85 
data.frame(
  率比 = unname(rate_result$estimate),
  CI下限 = unname(rate_result$conf.int[1]),
  CI上限 = unname(rate_result$conf.int[2]),
  p值 = rate_result$p.value,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "两样本精确 Poisson 率检验")
两样本精确 Poisson 率检验
率比 CI下限 CI上限 p值
新治疗 0.559 0.295 1.03 0.062

请核对 poisson.test() 的率比方向与输入顺序一致。简单 Poisson 方法假定事件按与人时相匹配的过程发生，且没有未处理的过度离散、患者内复发相关性或聚类。若同一患者可反复发生事件，或需调整年龄、中心和随访特征，应考虑 Poisson／负二项回归、稳健方差或专门的复发事件模型。

6 相关性检验

6.1 Pearson、Spearman 与 Kendall

相关分析描述两个变量共同变化的程度：

  • Pearson rr：衡量线性关联；对异常值敏感，区间和检验通常依赖成对观测独立及近似二元正态等条件。
  • Spearman ρ\rho：对数值取秩后衡量单调关联，适合有序资料或单调但非线性关系。
  • Kendall τ\tau：基于成对一致与不一致顺序，样本较小或并列秩较多时常有清晰解释。
plot(
  trial_data$age,
  trial_data$baseline_sbp,
  pch = 16,
  cex = 0.7,
  col = grDevices::adjustcolor(palette_test[["blue"]], alpha.f = 0.45),
  xlab = "年龄（岁）",
  ylab = "基线收缩压（mmHg）"
)
abline(
  lm(baseline_sbp ~ age, data = trial_data),
  col = palette_test[["vermillion"]],
  lwd = 2.5
)
散点图显示年龄与基线收缩压的关系，并叠加一条线性回归线。

基线年龄与收缩压的散点图及线性拟合。 

pearson_result <- cor.test(
  trial_data$age,
  trial_data$baseline_sbp,
  method = "pearson"
)
spearman_result <- cor.test(
  trial_data$age,
  trial_data$baseline_sbp,
  method = "spearman",
  exact = FALSE
)
kendall_result <- cor.test(
  trial_data$age,
  trial_data$baseline_sbp,
  method = "kendall",
  exact = FALSE
)

data.frame(
  方法 = c("Pearson", "Spearman", "Kendall"),
  相关估计 = c(
    unname(pearson_result$estimate),
    unname(spearman_result$estimate),
    unname(kendall_result$estimate)
  ),
  CI下限 = c(pearson_result$conf.int[1], NA, NA),
  CI上限 = c(pearson_result$conf.int[2], NA, NA),
  p值 = c(
    pearson_result$p.value,
    spearman_result$p.value,
    kendall_result$p.value
  ),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "三种常用相关分析")
三种常用相关分析
方法 相关估计 CI下限 CI上限 p值
Pearson 0.001 -0.109 0.111 0.985
Spearman -0.006 NA NA 0.918
Kendall -0.003 NA NA 0.935

base R 为 Pearson 相关直接给出常用置信区间；若需要秩相关区间，可使用与抽样设计匹配的自助法或专门软件。相关系数接近 0 只排除相应类型的关联，不排除所有关系。

相关不等于一致，也不等于因果 两台仪器即使 Pearson 相关接近 1，也可能存在固定偏差或随测量水平变化的偏差。评价连续测量一致性应结合差值图、Bland–Altman 界限或适当 ICC；分类判定一致性可考虑 kappa。相关还可能由混杂、范围限制、选择机制或共同时间趋势产生，不能单凭 p 值解释为因果关系。

7 生存结局：Kaplan–Meier 与 log-rank

7.1 为什么事件时间需要专门方法

生存资料同时包含事件是否发生、事件或删失发生的时间以及随时间变化的风险集。把删失者当作“无事件”，或只比较观察到事件者的平均时间，都会错误使用信息。

log-rank 检验比较两组完整生存曲线，在比例风险型差异下通常有较好效率。其零假设可理解为随访期间各事件时点的观察事件数与两组生存经历相同的预期相容。

survival_object <- survival::Surv(
  survival_data$time_months,
  survival_data$event
)
km_fit <- survival::survfit(
  survival_object ~ group,
  data = survival_data
)
logrank_result <- survival::survdiff(
  survival_object ~ group,
  data = survival_data
)
logrank_p <- pchisq(
  logrank_result$chisq,
  df = length(logrank_result$n) - 1,
  lower.tail = FALSE
)

data.frame(
  组别 = names(logrank_result$n),
  样本量 = unname(logrank_result$n),
  观察事件 = unname(logrank_result$obs),
  零假设期望事件 = unname(logrank_result$exp)
) |>
  knitr::kable(digits = 2, caption = "log-rank 检验的观察与期望事件数")
log-rank 检验的观察与期望事件数
组别 样本量.Var1 样本量.Freq 观察事件 零假设期望事件
group=标准治疗 A 130 86 63
group=新治疗 B 130 61 84 
data.frame(
  卡方统计量 = unname(logrank_result$chisq),
  自由度 = length(logrank_result$n) - 1,
  p值 = logrank_p,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "两组生存曲线的 log-rank 检验")
两组生存曲线的 log-rank 检验
卡方统计量 自由度 p值
15.1 1 0 
plot(
  km_fit,
  col = c(palette_test[["navy"]], palette_test[["vermillion"]]),
  lwd = 2.5,
  mark.time = TRUE,
  xlab = "随访时间（月）",
  ylab = "无事件概率",
  conf.int = FALSE
)
legend(
  "bottomleft",
  legend = levels(survival_data$group),
  col = c(palette_test[["navy"]], palette_test[["vermillion"]]),
  lwd = 2.5,
  bty = "n"
)
标准治疗与新治疗两组的 Kaplan-Meier 阶梯生存曲线，曲线上标有删失。

两组 Kaplan–Meier 生存曲线；短线表示删失。 

km_12 <- summary(km_fit, times = 12, extend = TRUE)
km_12_table <- data.frame(
  组别 = sub("group=", "", km_12$strata),
  时间月 = km_12$time,
  在险人数 = km_12$n.risk,
  无事件概率 = km_12$surv,
  CI下限 = km_12$lower,
  CI上限 = km_12$upper
)

knitr::kable(
  km_12_table,
  digits = 3,
  caption = "12 个月 Kaplan–Meier 无事件概率"
)
12 个月 Kaplan–Meier 无事件概率
组别 时间月 在险人数 无事件概率 CI下限 CI上限
标准治疗 12 40 0.423 0.342 0.522
新治疗 12 57 0.656 0.576 0.745 
cox_model <- survival::coxph(
  survival::Surv(time_months, event) ~ group,
  data = survival_data
)
cox_ci <- exp(confint(cox_model))
data.frame(
  比较 = "新治疗 vs 标准治疗",
  风险率比 = exp(coef(cox_model)),
  CI下限 = cox_ci[1, 1],
  CI上限 = cox_ci[1, 2],
  p值 = summary(cox_model)$coefficients[1, "Pr(>|z|)"],
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "未调整 Cox 模型的风险率比")
未调整 Cox 模型的风险率比
比较 风险率比 CI下限 CI上限 p值
group新治疗 新治疗 vs 标准治疗 0.523 0.375 0.73 0

log-rank 检验本身不提供效应大小，所以应同时给出 Kaplan–Meier 绝对生存概率、时间点、风险集和区间；若报告恒定 HR，则需检查比例风险假设。曲线明显交叉时，单个 log-rank p 值或 HR 可能难以概括差异，可考虑限制平均生存时间等与问题相符的估计量。

更多关于时间零点、删失、比例风险诊断、Cox PH 与 AFT 模型的内容，可继续阅读项目中的生存分析专题。

8 诊断试验与一致性

8.1 灵敏度、特异度与预测值

这里的“统计检验”不要与“诊断试验”混淆。评价一个检测需要预先定义目标人群、参考标准、阈值、时间关系以及不确定性。

diagnostic_table <- matrix(
  c(86, 14, 27, 173),
  nrow = 2,
  byrow = TRUE,
  dimnames = list(
    参考标准 = c("患病", "未患病"),
    新检测 = c("阳性", "阴性")
  )
)

tp <- diagnostic_table["患病", "阳性"]
fn <- diagnostic_table["患病", "阴性"]
fp <- diagnostic_table["未患病", "阳性"]
tn <- diagnostic_table["未患病", "阴性"]

sensitivity_ci <- binom.test(tp, tp + fn)$conf.int
specificity_ci <- binom.test(tn, tn + fp)$conf.int
ppv_ci <- binom.test(tp, tp + fp)$conf.int
npv_ci <- binom.test(tn, tn + fn)$conf.int

knitr::kable(
  addmargins(diagnostic_table),
  caption = "新检测与参考标准的交叉表"
)
新检测与参考标准的交叉表
阳性 阴性 Sum
患病 86 14 100
未患病 27 173 200
Sum 113 187 300 
diagnostic_metrics <- data.frame(
  指标 = c("灵敏度", "特异度", "阳性预测值", "阴性预测值"),
  估计值 = c(
    tp / (tp + fn),
    tn / (tn + fp),
    tp / (tp + fp),
    tn / (tn + fn)
  ),
  CI下限 = c(
    sensitivity_ci[1], specificity_ci[1], ppv_ci[1], npv_ci[1]
  ),
  CI上限 = c(
    sensitivity_ci[2], specificity_ci[2], ppv_ci[2], npv_ci[2]
  )
)

knitr::kable(
  diagnostic_metrics,
  digits = 3,
  caption = "诊断性能及精确二项 95% 置信区间"
)
诊断性能及精确二项 95% 置信区间
指标 估计值 CI下限 CI上限
灵敏度 0.860 0.776 0.921
特异度 0.865 0.810 0.909
阳性预测值 0.761 0.672 0.836
阴性预测值 0.925 0.878 0.958

灵敏度和特异度按参考标准分层；阳性和阴性预测值则以检测结果为条件，并随目标人群患病率变化。病例对照式抽样可估计灵敏度和特异度，但样本中的人为病例比例通常不能直接用于估计临床场景下的 PPV 和 NPV。

如果同一患者接受两种检测，比较阳性率、灵敏度或特异度时是配对设计。例如比较灵敏度时，应在参考标准确诊者中对两项检测的成对结果使用 McNemar 思路，而不是把两项检测当作两组独立样本。

9 检验前的数据检查与方法假设

9.1 先检查分母、编码、缺失与图形

最低限度的数据审核应包括：

  • 每组纳入人数、实际分析人数和缺失数；
  • 一位患者是否意外重复、配对标识是否唯一；
  • 单位、因子水平、事件编码和参考组是否正确；
  • 不可能值、录入错误、测量上限与下限；
  • 按组查看原始分布、异常值、散点图和时间轨迹；
  • 分析样本是否因完整案例筛选而改变目标人群。
data.frame(
  变量 = names(trial_data),
  类型 = vapply(trial_data, function(x) class(x)[1], character(1)),
  缺失数 = vapply(trial_data, function(x) sum(is.na(x)), integer(1)),
  唯一值数 = vapply(trial_data, function(x) length(unique(x)), integer(1)),
  check.names = FALSE
) |>
  knitr::kable(caption = "模拟试验数据的基本完整性检查")
模拟试验数据的基本完整性检查
变量 类型 缺失数 唯一值数
participant_id participant_id character 0 320
treatment treatment factor 0 2
age age numeric 0 56
sex sex factor 0 2
baseline_sbp baseline_sbp numeric 0 243
followup_sbp followup_sbp numeric 0 251
sbp_reduction sbp_reduction numeric 0 212
crp_week12 crp_week12 numeric 0 261
response response factor 0 2
adverse_event adverse_event factor 0 2 
stopifnot(!anyDuplicated(trial_data$participant_id))

异常值不能因为“删掉后显著”就删除。应先查明它是录入错误、测量错误、符合方案但罕见的真实值，还是代表目标人群之外的观察；保留或排除规则应基于数据质量和研究方案，并用敏感性分析展示影响。

9.2 正态性：图形与问题比单个 p 值更重要

t 检验或 ANOVA 的相关正态假设针对模型误差、组内分布或配对差值，而不是把所有组混在一起后的结局。应结合 Q–Q 图、样本量、偏态程度、异常值以及目标估计量判断。

anova_residuals <- residuals(traditional_anova)
qqnorm(
  anova_residuals,
  pch = 16,
  col = grDevices::adjustcolor(palette_test[["blue"]], alpha.f = 0.55),
  main = "ANOVA 残差 Q–Q 图"
)
qqline(anova_residuals, col = palette_test[["vermillion"]], lwd = 2.5)
正态分位数图展示三组方差分析残差与理论正态直线的接近程度。

三组 ANOVA 残差的正态 Q–Q 图。 

shapiro_result <- shapiro.test(anova_residuals)
fligner_result <- fligner.test(
  hemoglobin_change ~ dose,
  data = dose_data
)

data.frame(
  检查 = c("Shapiro–Wilk 残差正态性检验", "Fligner–Killeen 方差齐性检验"),
  统计量 = c(
    unname(shapiro_result$statistic),
    unname(fligner_result$statistic)
  ),
  p值 = c(shapiro_result$p.value, fligner_result$p.value),
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "两个常见假设检查检验")
两个常见假设检查检验
检查 统计量 p值
Shapiro–Wilk 残差正态性检验 0.99 0.235
Fligner–Killeen 方差齐性检验 9.05 0.011

不要把 Shapiro–Wilk 当作自动开关 小样本时该检验可能无法发现重要偏离，大样本时又可能对无实质影响的轻微偏离给出很小 p 值。更不应按“先检验正态；若不显著就 t 检验，否则 Wilcoxon”的两阶段规则自动选择，因为两个方法的目标估计量与假设并不相同。

同理，不必先用方差齐性检验决定是否允许 Welch t 检验；Welch 方法通常可直接作为两个独立均值比较的默认。Bartlett 检验对非正态敏感，Fligner–Killeen 更稳健，但任何方差检验都不能替代残差图和领域判断。

9.3 独立性通常无法从同一数据检验出来

独立性主要由抽样、随机化、随访和数据层级决定。若同一患者贡献两只眼、多个病灶或多次住院，仅看直方图无法发现错误的分析单位。应从数据字典与研究流程识别层级，并使用能表示相关结构的方法。

10 多重比较与选择性报告

10.1 为什么多次检验会累积假阳性风险

若在同一检验族中执行许多独立的 α=0.05\alpha=0.05 检验，至少出现一个假阳性的概率会增加。多重性不仅来自两两比较，还包括多个结局、时间点、亚组、阈值、模型版本和反复查看数据。

raw_p_values <- c(
  主要结局 = 0.012,
  次要结局A = 0.041,
  次要结局B = 0.007,
  次要结局C = 0.18,
  探索性结局 = 0.049
)

p_adjust_table <- data.frame(
  检验 = names(raw_p_values),
  原始p值 = unname(raw_p_values),
  Bonferroni = p.adjust(raw_p_values, method = "bonferroni"),
  Holm = p.adjust(raw_p_values, method = "holm"),
  BH_FDR = p.adjust(raw_p_values, method = "BH"),
  check.names = FALSE
)

knitr::kable(
  p_adjust_table,
  digits = 3,
  caption = "三种常见多重性调整示例"
)
三种常见多重性调整示例
检验 原始p值 Bonferroni Holm BH_FDR
主要结局 主要结局 0.012 0.060 0.048 0.030
次要结局A 次要结局A 0.041 0.205 0.123 0.061
次要结局B 次要结局B 0.007 0.035 0.035 0.030
次要结局C 次要结局C 0.180 0.900 0.180 0.180
探索性结局 探索性结局 0.049 0.245 0.123 0.061
  • Holm：控制检验族错误率，通常比单纯 Bonferroni 更有力，可作为一般确认性调整的实用选择。
  • Bonferroni：简单透明，但在检验较多或高度相关时可能保守。
  • Tukey：适合 ANOVA 中所有组间两两均值比较。
  • Dunnett 型比较：多组只与一个共同对照比较时更贴合目标。
  • Benjamini–Hochberg（BH）：控制错误发现率，常用于探索性或高维分析，不等同于控制至少一个假阳性的概率。

最有力的做法不是事后寻找最宽松的校正，而是在方案中明确主要结局、主要时间点、主要对比和同一检验族。调整 p 值时，置信区间也应与同一多重性策略一致。

11 优效性、等效性与非劣效性

11.1 未显著不等于相同

传统双侧优效性检验未拒绝“差异为零”，只表示数据对非零差异的证据不足。它不能证明两种治疗等效，也不能证明不存在临床重要差异。

等效性需要事先给定下、上等效界值，并证明效应区间完全落在该区间内。两个单侧检验（TOST）在显著性水平 α\alpha 下，等价于检查相应的 100(1−2α)%100(1-2\alpha)\% 置信区间是否完全位于等效界值之间。

mean_difference <- mean(new_values) - mean(standard_values)
standard_error <- sqrt(var(new_values) / n1 + var(standard_values) / n0)
welch_df <- (var(new_values) / n1 + var(standard_values) / n0)^2 /
  ((var(new_values) / n1)^2 / (n1 - 1) +
     (var(standard_values) / n0)^2 / (n0 - 1))

lower_equivalence_bound <- -6
upper_equivalence_bound <- 6
t_lower <- (mean_difference - lower_equivalence_bound) / standard_error
t_upper <- (mean_difference - upper_equivalence_bound) / standard_error
p_lower <- pt(t_lower, df = welch_df, lower.tail = FALSE)
p_upper <- pt(t_upper, df = welch_df, lower.tail = TRUE)
tost_p <- max(p_lower, p_upper)

ci90 <- mean_difference +
  c(-1, 1) * qt(0.95, df = welch_df) * standard_error

data.frame(
  均值差 = mean_difference,
  `90%CI下限` = ci90[1],
  `90%CI上限` = ci90[2],
  等效下界 = lower_equivalence_bound,
  等效上界 = upper_equivalence_bound,
  TOST_p值 = tost_p,
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "教学性 TOST 等效性示例")
教学性 TOST 等效性示例
均值差 90%CI下限 90%CI上限 等效下界 等效上界 TOST_p值
4.55 2.97 6.13 -6 6 0.066

等效界值必须来自临床判断、既往证据和方案，而不是为了让区间“刚好放进去”而在看完数据后设置。一个效应可能同时在统计上不同于 0，又在预先规定的宽等效区间内；“有差异”和“差异小到可视为等效”回答的是不同问题。

非劣效性则只关心一个预先规定方向：例如证明新治疗的疗效不比对照低超过界值 Δ\Delta。方向、界值、分析人群以及意向治疗和符合方案分析的角色都应事先定义。观察结果后改用单侧检验或重新选择界值会破坏错误率控制。

12 检验、回归模型与复杂设计

许多基础检验可写成回归模型：

基础检验 相应模型视角 扩展能力
两独立样本 t 检验 只有二分类组别的线性模型 加入基线值、协变量、交互和非线性
单因素 ANOVA 以组别因子为预测变量的线性模型 预设对比、不平衡设计、协变量调整
两比例／2×2 表分析 二元回归模型 调整混杂、估计 OR／RR／RD、交互
Poisson 率检验 带人时 offset 的 Poisson 模型 多变量率比、过度离散和聚类处理
log-rank 检验 未调整生存曲线比较 Cox、AFT、RMST 与协变量调整 
test_as_model <- lm(sbp_reduction ~ treatment, data = trial_data)
model_ci <- confint(test_as_model)["treatment新治疗", ]

data.frame(
  方法 = "线性模型：新治疗 vs 标准治疗",
  均值差 = coef(test_as_model)["treatment新治疗"],
  CI下限 = model_ci[1],
  CI上限 = model_ci[2],
  p值 = summary(test_as_model)$coefficients["treatment新治疗", "Pr(>|t|)"],
  check.names = FALSE
) |>
  knitr::kable(digits = 3, caption = "两组均值比较的线性模型表达")
两组均值比较的线性模型表达
方法 均值差 CI下限 CI上限 p值
treatment新治疗 线性模型：新治疗 vs 标准治疗 4.55 2.66 6.43 0

这里普通 lm() 的经典标准误对应等方差模型，因此数值不必与 Welch t 检验完全相同。通过合适的方差估计和模型结构可以放宽假设。

当研究包含多因素、中心效应、聚类、非线性、交互、删失或缺失随访时，单个基础检验往往不足。还要牢记：某亚组 p 值小、另一亚组 p 值大，并不证明效应在亚组间不同；应直接估计交互作用及其区间。

13 样本量、检验效能与临床意义

13.1 事前规划

样本量应围绕主要目标估计量、最小临床重要差异、预期变异、显著性水平、目标检验效能、分配比例、失访和设计效应事先规划。

continuous_power <- power.t.test(
  delta = 5,
  sd = 9,
  sig.level = 0.05,
  power = 0.80,
  type = "two.sample",
  alternative = "two.sided"
)

binary_power <- power.prop.test(
  p1 = 0.55,
  p2 = 0.70,
  sig.level = 0.05,
  power = 0.80,
  alternative = "two.sided"
)

data.frame(
  场景 = c("两独立均值差 5，SD 9", "两独立比例 0.55 vs 0.70"),
  每组所需完整资料人数 = ceiling(c(continuous_power$n, binary_power$n))
) |>
  knitr::kable(caption = "简化假设下的样本量示例")
简化假设下的样本量示例
场景 每组所需完整资料人数
两独立均值差 5，SD 9 52
两独立比例 0.55 vs 0.70 163

这些函数适用于简化的独立两组设计。配对、聚类、重复测量、生存结局、不等分配、多个主要结局或复杂模型需要专门公式或仿真。计算出的完整资料人数还应按合理失访率膨胀。

不推荐在结果出来后用观察到的效应计算“事后功效”来解释不显著结果，因为它通常只是 p 值的重新表达。此时应查看效应估计和 95% CI：区间是否排除了临床重要获益、伤害或两者都没有？

14 正确理解 p 值与置信区间

14.1 p 值能说什么，不能说什么

p 值是在零假设及统计模型假设成立时，得到当前观察结果或更不利于零假设结果的概率。它不是：

  • 零假设为真的概率；
  • 结果“由偶然造成”的概率；
  • 效应大小或临床重要性的度量；
  • 研究没有偏倚、混杂或测量问题的证明；
  • 另一项研究必然重复成功的概率。

固定 α=0.05\alpha=0.05 时，第一类错误是零假设成立却被拒绝；第二类错误是有研究所定义的效应却未拒绝零假设。检验效能是给定真实效应和设计条件下正确拒绝零假设的概率。

置信区间把效应大小和不确定性放在同一尺度。95% 置信区间的频率学解释是：若按相同设计和程序反复抽样，长期来看约 95% 的这样构造的区间会覆盖真实参数。它不是“这个已得到的区间有 95% 概率包含真实值”。

14.2 单侧与双侧检验

医学研究通常使用双侧检验，因为未预期方向的伤害或反向效应同样重要。只有当反方向效应在科学和决策上确实无需区分、方向在看数据前已写入方案且监管或领域规范允许时，单侧检验才可能合理。观察结果后把双侧改成单侧相当于选择性降低 p 值。

更好的阅读顺序：先看每组分母与描述统计，再看预先规定方向的效应估计及 95% CI，随后看 p 值和多重性调整，最后判断设计、数据质量、临床意义与可推广性。

15 如何报告检验结果

15.1 最低报告清单

一段可审核的医学统计结果至少应说明：

  • 研究设计、分析单位、组别或配对方式；
  • 每组实际分析人数、缺失数和关键分母；
  • 与结局尺度匹配的描述统计；
  • 对比方向、效应量、单位和 95% CI；
  • 检验或模型的完整名称，单侧或双侧；
  • 精确 p 值（极小时写 p < 0.001，不要写 p = 0）；
  • 多重性处理、假设检查和任何敏感性分析；
  • 结果适用范围以及设计、偏倚或模型限制。

15.2 可复用报告模板

在【人群与设计】中，【组 A】与【组 B】的【结局】分别为【带分母的描述统计】。预先指定方向的【效应尺度】为【估计值与单位】，95% CI【下限，上限】；使用【完整检验名称、单双侧及校正方法】得到 p=p=【数值】。分析基于【实际样本量】份完整／可用记录，并受【具体假设、缺失、偏倚或推广限制】影响。

15.3 本教程连续结局示例

report_difference <- mean(new_values) - mean(standard_values)
report_ci <- -rev(welch_result$conf.int)
report_p <- format_p(welch_result$p.value)

cat(
  paste0(
    "新治疗组与标准治疗组各纳入 ", n1, " 与 ", n0,
    " 人，平均收缩压下降分别为 ",
    sprintf("%.1f", mean(new_values)), "（SD ",
    sprintf("%.1f", sd(new_values)), "）和 ",
    sprintf("%.1f", mean(standard_values)), "（SD ",
    sprintf("%.1f", sd(standard_values)), "）mmHg。",
    "新治疗减标准治疗的均值差为 ",
    sprintf("%.1f", report_difference), " mmHg（95% CI ",
    sprintf("%.1f", report_ci[1]), " 至 ",
    sprintf("%.1f", report_ci[2]), "；双侧 Welch t 检验 p ",
    ifelse(welch_result$p.value < 0.001, "< 0.001", paste0("= ", report_p)),
    "）。该估计来自模拟完整资料，临床解释还需结合预先规定的最小临床重要差异。"
  )
)
## 新治疗组与标准治疗组各纳入 160 与 160 人，平均收缩压下降分别为 11.6（SD 8.9）和 7.1（SD 8.2）mmHg。新治疗减标准治疗的均值差为 4.5 mmHg（95% CI 2.7 至 6.4；双侧 Welch t 检验 p < 0.001）。该估计来自模拟完整资料，临床解释还需结合预先规定的最小临床重要差异。

16 常见错误速查

错误做法 为什么有问题 更合适的处理
p > 0.05 就写“两组相同” 未拒绝零假设不证明相同 报告效应与 CI；等效问题用预设界值
只写“差异有统计学意义” 没有方向、大小、单位和精度 报告原始尺度效应、CI 与临床意义
先做 Shapiro 检验再机械选方法 两阶段规则忽略目标估计量与稳健性 先定问题和设计，结合图形及敏感性分析
把 Mann–Whitney 称为普遍的中位数检验 一般检验的是秩分布 说明额外形状假设或使用准确表述
把配对数据当独立样本 忽略患者内相关性 使用配对检验或重复测量模型
对三个时点执行三次未校正检验 增加检验族假阳性率 先总体模型，再做预设且校正的对比
因异常值影响显著性便删除 形成结果驱动的数据处理 查明来源，按方案处理并做敏感性分析
把 OR 当 RR 结局常见时二者可差很多 按设计报告 RD、RR 或 OR，并写清尺度
RCT 基线变量逐项做显著性检验 随机化后差异本来就是随机产生 描述基线；按方案基于预后价值调整
一个亚组显著、另一个不显著就称有交互 两个 p 值的差不等于差异的 p 值 直接估计交互项及其 CI
看结果后改成单侧检验 破坏预设错误率 在方案中事先决定方向和理由
隐藏不同分析的样本量 完整案例可能改变样本与目标人群 每项分析报告实际分母和缺失
不断更换检验、阈值或结局直到显著 产生选择性报告与不稳定结论 预注册主要分析并完整披露探索过程
把相关当成一致性或因果 三者回答不同问题 使用一致性方法，并基于设计讨论因果

17 练习与答案

17.1 练习 1：两个独立均值

一项随机试验比较两组 8 周 LDL-C 变化，每位参与者只属于一组。两组方差看起来不完全相同，研究问题是平均变化差。最直接的基础检验是什么？主要报告什么？

查看答案 使用 Welch 两独立样本 t 检验。主要报告每组样本量、均值和 SD，以及预先规定方向的平均变化差、95% CI 和双侧 p 值。不能因为 Welch p 值显著就省略临床单位下的差异。

17.2 练习 2：同一患者的左右眼

研究者比较同一患者左眼和右眼的眼压，却使用独立样本 t 检验。核心错误是什么？

查看答案 左右眼属于同一患者，观测不是独立样本。若每位患者恰好有一对眼值且问题是左右差异，可分析患者内差值；更复杂的眼别、治疗、重复访视或双眼入组设计应使用能够表示眼嵌套于患者相关性的模型。

17.3 练习 3：三组总体检验

三组 ANOVA 的 p 值为 0.003。能否直接写“三组彼此都不同”？

查看答案 不能。总体检验只表明至少有一组均值不符合全部相等的零假设。应估计预先指定的对比，或使用 Tukey、Dunnett、Holm 等与比较族相匹配的方法，并报告同时区间或调整后的推断。

17.4 练习 4：稀疏不良事件

两组各 25 人，一组 0 例严重不良事件，另一组 3 例。为什么不能只报告 Fisher 检验 p 值？

查看答案 小样本和零事件使效应估计非常不精确。还应报告每组分母与风险、风险差或适当相对效应及区间，并说明零单元导致的估计问题。未显著不能排除临床重要伤害。

17.5 练习 5：前后阳性率

同一批患者在干预前后接受同一二分类检测。哪个检验利用了正确的配对结构？

查看答案 McNemar 检验；若不一致配对很少，可对两个方向的不一致对使用精确二项检验。还应报告前后边际比例及风险差，而不是只给 p 值。

17.6 练习 6：相关性很高

新仪器与标准仪器的 Pearson r=0.97r=0.97。这是否足以证明两台仪器可以互换？

查看答案 不足。高度相关可以在存在固定或比例偏差时出现。应检查配对差值随测量水平的模式，报告 Bland–Altman 偏差与一致性界限，并根据用途考虑 ICC、重复性和可接受误差范围。

17.7 练习 7：生存曲线交叉

两组 Kaplan–Meier 曲线在 6 个月交叉，之后差异方向反转。为什么单个 log-rank p 值可能不够？

查看答案 交叉提示效应随时间变化，log-rank 在比例风险型差异下的优势和单个恒定 HR 的解释都可能受限。应展示曲线、风险表和关键时点绝对生存概率，检查时间变化效应，并考虑与研究问题相符的 RMST 或分段效应。

17.8 练习 8：不显著结果

治疗效应估计为风险差 −2%-2\%，95% CI 为 −18%-18\% 至 14%14\%，p=0.81p=0.81。能否得出“没有治疗作用”？

查看答案 不能。区间同时包含可能有意义的获益和伤害，说明估计不精确。可说数据未显示明确差异，但不能证明无作用；应讨论样本量、事件数、临床重要界值和研究偏倚。

18 快速参考

18.1 常用 R 代码

分析目的 base R 或推荐包代码模式 优先报告
单样本均值 t.test(x, mu = value) 均值与参考值之差、95% CI
两独立均值 t.test(y ~ group) 均值差、95% CI
配对均值 t.test(after, before, paired = TRUE) 平均个体差、95% CI
两独立秩分布 wilcox.test(y ~ group, conf.int = TRUE) 分布、位置差与 CI
配对秩分布 wilcox.test(after, before, paired = TRUE) 差值分布、位置差
多组均值 oneway.test(y ~ group, var.equal = FALSE) 预设组间均值差
传统 ANOVA aov(y ~ group) 总体检验、具体对比与效应量
多组秩分布 kruskal.test(y ~ group) 总体秩差与校正后比较
完整区组重复测量 friedman.test(y ~ time \| id) 时点差异与校正后比较
单比例 binom.test(x, n) 比例与精确 CI
两独立比例 prop.test(x, n)；chisq.test(tab) RD、RR 或 OR 与 CI
稀疏 2×2 表 fisher.test(tab) 分母、绝对风险、OR 与 CI
配对二分类 mcnemar.test(tab) 边际风险差、不一致对
有序比例趋势 prop.trend.test(x, n) 各级比例与趋势检验
分层 2×2 表 mantelhaen.test(array) 分层结果、共同 OR 与 CI
两个发病率 poisson.test(events, person_time) 率、率比与 CI
Pearson 相关 cor.test(x, y, method = "pearson") rr 与 CI
秩相关 cor.test(x, y, method = "spearman") ρ\rho 及适当区间
生存曲线比较 survival::survdiff(Surv(time, event) ~ group) KM 概率、风险集、效应与 CI
多重 p 值调整 p.adjust(p, method = "holm") 调整方法、调整 p 与同时 CI

18.2 零效应值速查

效应尺度 零效应值 典型解释
均值差、风险差、率差、相关系数 0 无相应尺度上的差异或关联
风险比、率比、优势比、风险率比 1 两组相对尺度相同
灵敏度、特异度、预测值 无统一“零效应”值 与用途、阈值和基准比较

最终分析检查清单

发布结果前，请确认：

  • 研究问题、总体、时间窗和目标估计量已预先写清；
  • 分析单位以及独立、配对、重复、聚类或删失结构识别正确；
  • 结局编码、对比方向、单位、参考组和分母均已核对；
  • 数据质量、缺失、异常值和各组分布已经图形化检查；
  • 方法假设与研究设计相符，且没有由单个正态性 p 值机械决定方法；
  • 每项主要结论都有原始尺度效应量与 95% CI；
  • p 值的单双侧方向、精确方法和调整方式已说明；
  • 多结局、多时点、多亚组和模型选择产生的多重性已处理或透明披露；
  • 统计显著没有被误写成临床重要，未显著没有被误写成相同；
  • OR、RR、HR、率比、相关与因果等概念没有混用；
  • 实际分析样本量、缺失记录和敏感性分析均已报告；
  • 代码、随机种子、软件版本和数据处理轨迹足以支持复现。
sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/Edmonton
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   Matrix_1.7-5   
##  [5] xfun_0.60       lattice_0.22-9  splines_4.6.1   cachem_1.1.0   
##  [9] knitr_1.51      htmltools_0.5.9 rmarkdown_2.31  lifecycle_1.0.5
## [13] cli_3.6.6       grid_4.6.1      sass_0.4.10     jquerylib_0.1.4
## [17] compiler_4.6.1  tools_4.6.1     evaluate_1.0.5  bslib_0.12.0
## [21] survival_3.8-6  yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0

本教程中的数据与结果全部为模拟内容，仅用于医学统计学教学，不能替代研究方案阶段的统计咨询、临床判断、伦理审查或监管要求。真实研究应根据设计、抽样、数据质量和目标估计量选择分析方法。