The V Lab
适用对象医学、心理学、药学与公共卫生学习者
学习时长约 180–240 分钟
先修要求置信区间、回归模型与基础 R

教学案例与责任边界 本页的 360 名参与者、研究中心、疗效和不良事件均由固定随机种子模拟,不对应任何真实患者或产品。代码用来解释设计与分析原则,不能代替试验方案、统计分析计划、伦理审查、数据监查委员会或经过验证的监管分析流程。

如何使用本教程

临床试验不是“收完数据再选检验”。建议按下面的因果与操作链条学习:

临床问题 → estimand → 设计保护 → 样本量 → 数据质量 → 预设模型 → 敏感性分析 → 透明报告

贯穿案例是一项模拟的多中心、双盲、两臂平行优效性试验。主结局是 12 周症状评分,分数越低越好;同时演示应答、急性发作次数、不良事件和一年内复发。所有效应统一写成 Treatment − Control:连续结局负值有利,应答的 RD/RR 大于基准有利,不良事件的 RR 大于 1 表示风险增加,复发 HR 小于 1 有利。

学习目标

完成本教程后,你应能够:

  • 区分试验阶段、研究目的和具体设计;
  • 把 PICO 问题写成 ICH E9(R1) 的五要素 estimand;
  • 区分随机序列、分配隐藏和盲法,评价基线平衡而不“检验随机化”;
  • 从主要 estimand 倒推连续、二元和时间结局的样本量;
  • 为连续、重复测量、二元、计数和生存结局选择模型并报告效应量与 95% CI;
  • 区分 ITT/FAS、符合方案、安全性和实际治疗分析集;
  • 处理随机化分层因素、协变量调整、缺失数据与 intercurrent events;
  • 正确解释多重性、期中分析、亚组交互和非劣效界值;
  • 为平行、整群、交叉、析因和适应性试验识别关键分析单位;
  • 按 CONSORT 2025 报告参与者流程、分析集、效应、不确定性与伤害。

1 临床试验回答什么问题

1.1 介入性问题与观察性问题

临床试验由研究者分配干预,目标通常是估计一个明确治疗策略的因果效应;观察研究记录已发生的暴露,需要额外处理混杂。随机化的核心价值是让已知和未知基线预后因素在概率上可比,而不是保证每个小样本的每一项变量数值相等。

“解释性”试验更关注理想条件下能否起效(efficacy),“务实性”试验更接近日常路径中的效果(effectiveness)。两者是连续谱:入选标准、干预弹性、对照、随访强度和分析策略都应服务于目标使用场景。

1.2 阶段不等于统计模型

阶段 主要目标 常见重点 统计提醒
Phase 0 / 早期探索 微剂量、靶点或可行性 PK/PD、机制 通常不用于确证疗效
Phase I 安全性、耐受性、剂量 DLT、PK/PD、剂量递增 小样本不代表“无需设计”
Phase II 剂量选择与初步活性 剂量–反应、可行终点 选择规则会影响后续不确定性
Phase III 确证获益–风险 预设主要结局、对照、多中心 控制 I 类错误并保护估计量
Phase IV 上市后有效性与安全性 罕见伤害、真实世界实施 长期、依从性与选择问题突出

器械、行为干预、心理治疗和罕见病研究未必沿用同一种药物阶段标签。阶段说明开发目标,设计与 estimand 才决定分析。

1.3 方案、注册与统计分析计划

在看到分组结果以前,至少应锁定:入排标准、时间零点、干预版本、对照、主要和关键次要结局、测量时间、estimand、随机化、样本量、分析集、协变量、缺失数据、多重性、期中规则和敏感性分析。

  • Protocol:说明为什么和如何开展试验;可按 SPIRIT 2025 组织。
  • 注册记录:让研究问题、主要结局和时间线可追踪;注册不等于完整方案。
  • SAP:在揭盲或访问非盲分组结果前细化变量推导、模型、对比、异常情况和输出。
  • 修订:允许有理由的变更,但必须带日期、理由、责任人,并区分揭盲前后。

2 从 PICO 到 estimand

2.1 五个不可缺的属性

ICH E9(R1) 要求把治疗问题写成 estimand。贯穿案例的一个可操作版本如下:

属性 本案例的主要 estimand
Population 符合条件且被随机分配的成人,多中心目标人群
Treatment conditions 分配 Treatment 相对分配 Control,包括方案规定的伴随照护
Variable 从基线至第 12 周的症状评分;第 12 周分数为主要分析变量
Intercurrent-event strategy 停药或换救援治疗采用 treatment-policy;死亡等无法定义量表的事件需另定策略
Population-level summary 调整基线及随机化分层因素后的第 12 周平均差(Treatment − Control)

常见 intercurrent-event 策略不是缺失值“补法”的同义词:

  • Treatment policy:不论停药、换药等事件,比较原始分配策略的结果;
  • Hypothetical:估计如果特定事件没有发生会怎样;
  • Composite:把事件本身并入结局,例如停药即判无应答;
  • While on treatment:只关心事件发生前的结局;
  • Principal stratum:针对在两种潜在治疗下都满足某事件状态的人群,通常需要强假设。

ICE 与缺失不是一回事 停药、救援用药、死亡或妊娠是随机化后的事件;量表未测到是缺失数据。一个 intercurrent event 可以导致缺失,但必须先定义临床问题,再选择与该问题一致的数据处理和估计方法。

2.2 主分析、敏感性与补充分析

估计量(estimator)是从数据得到 estimand 估计值的规则。主分析依赖一组明确假设;敏感性分析改变关键且不可验证的假设(例如把治疗组缺失值系统性调差);补充分析回答相邻但不同的问题(例如符合方案效果)。“多跑几个模型看是否显著”不属于预设敏感性框架。

3 设计保护:随机化、隐藏与盲法

3.1 三件事分别解决三种偏倚

要素 问题 实务要求
随机序列生成 下一位应分到哪一组? 可复现的算法、适当分层/区组、受控随机化系统
分配隐藏 入组前能否预测或改变分组? 招募者不能看到后续序列,变区组大小并限制访问
盲法 分组后行为或评价是否受知晓影响? 分别说明参与者、治疗者、结局评价者和分析者

盲法并非一句“双盲”即可。某些手术或心理干预无法让治疗者盲,但仍可采用独立盲态结局评价、标准化共同干预、盲态终点判定和预设分析。紧急揭盲的触发、权限和记录也属于设计。

3.2 分层变区组随机化的模拟

下面生成 360 名模拟参与者。以研究中心 × 基线严重程度分层,在每层使用大小为 4 或 6 的随机区组;区组大小在真实试验中不应向招募人员公开。

set.seed(20260811)

make_block_sequence <- function(n, block_sizes = c(4L, 6L)) {
  allocation <- character(0)
  while (length(allocation) < n) {
    block_size <- sample(block_sizes, 1L)
    block <- rep(c("Control", "Treatment"), each = block_size / 2L)
    allocation <- c(allocation, sample(block))
  }
  allocation[seq_len(n)]
}

n <- 360L
trial <- data.frame(
  id = sprintf("P%03d", seq_len(n)),
  site = factor(sample(paste0("Site ", 1:4), n, replace = TRUE,
                       prob = c(0.28, 0.26, 0.24, 0.22))),
  sex = factor(sample(c("Female", "Male"), n, replace = TRUE,
                      prob = c(0.58, 0.42))),
  age = round(pmin(pmax(rnorm(n, 45, 13), 18), 75), 1)
)

site_baseline <- c("Site 1" = 0, "Site 2" = 0.8,
                   "Site 3" = -0.6, "Site 4" = 0.4)
trial$baseline <- round(
  pmin(pmax(rnorm(n, 24, 5.5) + site_baseline[trial$site], 10), 40), 1
)
trial$severity <- factor(ifelse(trial$baseline >= 24, "High", "Low"),
                         levels = c("Low", "High"))
trial$stratum <- interaction(trial$site, trial$severity, drop = TRUE)

trial$arm <- NA_character_
for (stratum_i in levels(trial$stratum)) {
  index <- which(trial$stratum == stratum_i)
  trial$arm[index] <- make_block_sequence(length(index))
}
trial$arm <- factor(trial$arm, levels = c("Control", "Treatment"))
treatment <- as.integer(trial$arm == "Treatment")

  # 第 12 周连续结局:分数越低越好。
site_week12 <- c("Site 1" = 0, "Site 2" = 0.5,
                 "Site 3" = -0.5, "Site 4" = 0.2)
trial$week12_full <- 5 + 0.58 * trial$baseline - 2.5 * treatment +
  site_week12[trial$site] + 0.25 * (trial$sex == "Male") + rnorm(n, 0, 5.2)
trial$week12_full <- round(pmin(pmax(trial$week12_full, 0), 45), 1)

  # 较差结局和治疗组有更高的缺失概率;完整值只用于教学模拟真值。
p_missing <- plogis(-2.25 + 0.35 * treatment +
                     0.085 * (trial$week12_full - 18) +
                     0.20 * (trial$severity == "High"))
trial$missing12 <- rbinom(n, 1, p_missing)
trial$week12 <- ifelse(trial$missing12 == 1, NA_real_, trial$week12_full)
trial$change <- trial$week12 - trial$baseline

  # 应答:较基线下降至少 30%;缺失按无应答处理只是一个 composite 策略示例。
trial$responder_full <- as.integer(
  (trial$baseline - trial$week12_full) / trial$baseline >= 0.30
)
trial$responder_nri <- ifelse(is.na(trial$week12), 0L, trial$responder_full)

  # 安全性结局。
p_ae <- plogis(qlogis(0.13) + log(1.75) * treatment +
                0.012 * (trial$age - 45))
trial$adverse_event <- rbinom(n, 1, p_ae)

  # 一年内复发/进展与独立失访删失。
event_time <- rexp(n, rate = 0.0022 *
                    exp(log(0.68) * treatment + 0.020 * (trial$baseline - 24)))
dropout_time <- rexp(n, rate = 0.00075)
trial$time_days <- pmin(event_time, dropout_time, 365)
trial$status <- as.integer(event_time <= pmin(dropout_time, 365))

  # 重复评分与带过度离散的 12 周发作计数;放在核心结局后模拟以固定主结果。
trial$week4 <- pmax(0, 3.5 + 0.78 * trial$baseline - 0.9 * treatment +
                    rnorm(n, 0, 4.7))
trial$week8 <- pmax(0, 4.0 + 0.66 * trial$baseline - 1.8 * treatment +
                    rnorm(n, 0, 5.0))
trial$week4[rbinom(n, 1, plogis(-3 + 0.2 * treatment)) == 1] <- NA
trial$week8[rbinom(n, 1, plogis(-2.7 + 0.25 * treatment)) == 1] <- NA
trial$exposure_weeks <- round(runif(n, 8, 12), 1)
individual_frailty <- rgamma(n, shape = 2, rate = 2)
episode_mean <- exp(log(0.16) + log(0.75) * treatment +
                    0.018 * (trial$baseline - 24)) *
  trial$exposure_weeks * individual_frailty
trial$episode_count <- rpois(n, episode_mean)

stopifnot(nrow(trial) == n, !anyNA(trial$arm),
          all(trial$status %in% 0:1), all(trial$time_days > 0))

allocation_table <- addmargins(table(trial$stratum, trial$arm))
allocation_table
##              
##               Control Treatment Sum
##   Site 1.Low       22        24  46
##   Site 2.Low       18        19  37
##   Site 3.Low       28        28  56
##   Site 4.Low       16        16  32
##   Site 1.High      27        28  55
##   Site 2.High      28        27  55
##   Site 3.High      17        18  35
##   Site 4.High      22        22  44
##   Sum             178       182 360

3.3 基线表用描述与 SMD,不做显著性狩猎

随机化后基线差异是随机变量。逐项 p 值不能证明随机化成功,也不应决定是否调整;调整变量应由预后价值、分层设计和 SAP 预设。标准化差异(SMD)帮助看量级,但也不是“通过/不通过”检验。

smd_continuous <- function(x, arm) {
  x0 <- x[arm == "Control"]
  x1 <- x[arm == "Treatment"]
  (mean(x1) - mean(x0)) /
    sqrt((stats::var(x1) + stats::var(x0)) / 2)
}

smd_binary <- function(x, arm) {
  p0 <- mean(x[arm == "Control"])
  p1 <- mean(x[arm == "Treatment"])
  (p1 - p0) / sqrt((p1 * (1 - p1) + p0 * (1 - p0)) / 2)
}

baseline_table <- data.frame(
  variable = c("样本量", "基线评分,均值 (SD)", "年龄,均值 (SD)",
               "男性,n (%)", "第 12 周缺失,n (%)"),
  Control = c(
    sum(trial$arm == "Control"),
    sprintf("%.1f (%.1f)", mean(trial$baseline[trial$arm == "Control"]),
            sd(trial$baseline[trial$arm == "Control"])),
    sprintf("%.1f (%.1f)", mean(trial$age[trial$arm == "Control"]),
            sd(trial$age[trial$arm == "Control"])),
    sprintf("%d (%.1f%%)", sum(trial$sex == "Male" & trial$arm == "Control"),
            100 * mean(trial$sex[trial$arm == "Control"] == "Male")),
    sprintf("%d (%.1f%%)", sum(trial$missing12 == 1 & trial$arm == "Control"),
            100 * mean(trial$missing12[trial$arm == "Control"]))
  ),
  Treatment = c(
    sum(trial$arm == "Treatment"),
    sprintf("%.1f (%.1f)", mean(trial$baseline[trial$arm == "Treatment"]),
            sd(trial$baseline[trial$arm == "Treatment"])),
    sprintf("%.1f (%.1f)", mean(trial$age[trial$arm == "Treatment"]),
            sd(trial$age[trial$arm == "Treatment"])),
    sprintf("%d (%.1f%%)", sum(trial$sex == "Male" & trial$arm == "Treatment"),
            100 * mean(trial$sex[trial$arm == "Treatment"] == "Male")),
    sprintf("%d (%.1f%%)", sum(trial$missing12 == 1 & trial$arm == "Treatment"),
            100 * mean(trial$missing12[trial$arm == "Treatment"]))
  ), check.names = FALSE
)

smd_table <- data.frame(
  variable = c("基线评分", "年龄", "男性"),
  SMD = c(smd_continuous(trial$baseline, trial$arm),
          smd_continuous(trial$age, trial$arm),
          smd_binary(trial$sex == "Male", trial$arm))
)

knitr::kable(baseline_table, caption = "模拟试验的基线与随访完整性")
模拟试验的基线与随访完整性
variable Control Treatment
样本量 178 182
基线评分,均值 (SD) 24.2 (5.4) 24.4 (6.0)
年龄,均值 (SD) 45.1 (13.1) 44.9 (13.2)
男性,n (%) 78 (43.8%) 73 (40.1%)
第 12 周缺失,n (%) 20 (11.2%) 25 (13.7%)
knitr::kable(smd_table, digits = 3,
             caption = "Treatment − Control 标准化差异;绝对值仅描述量级")
Treatment − Control 标准化差异;绝对值仅描述量级
variable SMD
基线评分 0.030
年龄 -0.016
男性 -0.075

区组随机化使总人数接近 1:1,三个 SMD 的绝对值均很小。但这一观察结果不是随机化质量的主要证据;更关键的是序列生成、隐藏、实施日志和分配后排除是否可审计。

4 优效、非劣效与等效

4.1 先确定科学框架,再看置信区间

设连续结局越低越好,效应为 θ=μT−μC\theta=\mu_T-\mu_C:

目标 原假设的核心 支持结论的 95% CI 位置 不能使用的捷径
优效性 治疗不优于对照 双侧 CI 排除 0,方向有利 只说“治疗组内显著”
非劣效性 治疗比对照差至少 Δ\Delta 单侧 95% 上界低于 +Δ+\Delta “差异不显著,所以不劣”
等效性 差异在容许区间之外 双侧 90% CI 完全落在 [−Δ,+Δ][-\Delta,+\Delta] 把高 p 值当等效证据

非劣效界值必须在看结果前由临床可接受损失与可靠历史证据共同论证,还要讨论活性对照在当前试验能起效的 assay sensitivity、历史效应可迁移的 constancy 假设和连续“放宽界值”导致的 biocreep。非劣效试验通常并列 FAS/ITT 与 PP 结果;一致性增强解释,但不能自动消除两者共同的偏倚。

  # 教学演示:沿用优效性模拟数据来说明 CI 判定,不构成真实 NI 设计。
fit_ancova_preview <- lm(
  week12 ~ arm + baseline + site + severity,
  data = trial
)
ni_margin <- 2
ni_estimate <- coef(fit_ancova_preview)["armTreatment"]
ni_se <- sqrt(vcov(fit_ancova_preview)["armTreatment", "armTreatment"])
ni_upper <- ni_estimate + qt(0.95, df.residual(fit_ancova_preview)) * ni_se

ni_table <- data.frame(
  comparison = "Treatment − Control(低分有利)",
  estimate = ni_estimate,
  one_sided_95_upper = ni_upper,
  margin = ni_margin,
  noninferior = ni_upper < ni_margin
)
knitr::kable(ni_table, digits = 3,
             caption = "非劣效 CI 逻辑演示;界值 +2 仅为教学设定")
非劣效 CI 逻辑演示;界值 +2 仅为教学设定
comparison estimate one_sided_95_upper margin noninferior
armTreatment Treatment − Control(低分有利) -1.76 -0.775 2 TRUE

5 常见设计及其分析单位

设计 适用场景 关键分析 常见陷阱
个体平行组 大多数急慢性问题 按随机分组比较,调整预设分层因素 分配后排除、选择性换结局
整群随机 污染明显或干预在机构层实施 用整群相关结构;样本量含 ICC 与整群数 把个体当独立、整群太少
交叉 病情稳定、作用可逆且洗脱可行 纳入 period、sequence、个体内相关 残留效应、不可逆结局
析因 同时研究两个干预 主效应与交互按 estimand 预设 有交互却只解释主效应
适应性/平台 预设的剂量、样本量或臂调整 模拟 operating characteristics,控制 I 类错误 临时规则、操作偏倚、时间趋势

整群设计的简单设计效应为 DE=1+(m‾−1)ICCDE=1+(\bar m-1)ICC,但整群大小不等、分层匹配和少量整群会使公式不足。交叉、析因和适应性设计的更多推导见实验设计专题;试验中必须另外保护分配隐藏、患者安全、可解释 estimand 和监管审计轨迹。

适应性不等于边做边改 适应规则、信息时间、决策者、模拟、错误率与数据访问必须预设。当前可参考适应性设计原则和 draft ICH E20;不要把尚未最终采纳的草案写成已生效标准。

6 样本量:从主要 estimand 倒推

6.1 样本量不是一个软件按钮

计算前应明确主要结局的尺度、目标差异、方差或对照风险、双侧/单侧 α\alpha、把握度、分配比、失访、协变量效率、整群或重复测量相关、期中设计和多重性。目标差异应有临床意义,不能只取预试验中最乐观的点估计。

连续结局、两组等分配的近似为:

n每组≈2σ2(z1−α/2+z1−β)2δ2. n_{\text{每组}}\approx \frac{2\sigma^2(z_{1-\alpha/2}+z_{1-\beta})^2}{\delta^2}.

二元结局需要两组风险而不只是相对效应;时间到事件试验主要由事件数驱动;整群试验还取决于整群数与 ICC。失访膨胀通常用 n/(1−d)n/(1-d),不是 n(1+d)n(1+d)。

continuous_ss <- power.t.test(
  delta = 2.5, sd = 6, sig.level = 0.05, power = 0.80,
  type = "two.sample", alternative = "two.sided"
)
binary_ss <- power.prop.test(
  p1 = 0.35, p2 = 0.50, sig.level = 0.05, power = 0.80,
  alternative = "two.sided"
)

sample_size_table <- data.frame(
  endpoint = c("连续评分:差 2.5,SD 6", "应答:35% 对 50%"),
  per_arm = c(ceiling(continuous_ss$n), ceiling(binary_ss$n)),
  per_arm_with_15pct_attrition = c(
    ceiling(continuous_ss$n / 0.85), ceiling(binary_ss$n / 0.85)
  )
)
knitr::kable(sample_size_table,
             caption = "80% 把握度、双侧 α=0.05、1:1 分配的教学计算")
80% 把握度、双侧 α=0.05、1:1 分配的教学计算
endpoint per_arm per_arm_with_15pct_attrition
连续评分:差 2.5,SD 6 92 108
应答:35% 对 50% 170 200
scenario_grid <- expand.grid(
  target_difference = c(2.0, 2.5, 3.0),
  sd = c(5, 6),
  attrition = c(0.10, 0.20)
)
scenario_grid$per_arm <- mapply(function(delta, sd, loss) {
  base_n <- power.t.test(delta = delta, sd = sd, sig.level = 0.05,
                         power = 0.80, type = "two.sample")$n
  ceiling(base_n / (1 - loss))
}, scenario_grid$target_difference, scenario_grid$sd, scenario_grid$attrition)
knitr::kable(scenario_grid, caption = "关键假设变化时的连续结局样本量情景")
关键假设变化时的连续结局样本量情景
target_difference sd attrition per_arm
2.0 5 0.1 111
2.5 5 0.1 71
3.0 5 0.1 50
2.0 6 0.1 159
2.5 6 0.1 102
3.0 6 0.1 71
2.0 5 0.2 124
2.5 5 0.2 80
3.0 5 0.2 56
2.0 6 0.2 178
2.5 6 0.2 115
3.0 6 0.2 80

不要报告事后把握度 观察到数据后,“observed power”基本是 p 值的重新编码,不能补充效应估计和置信区间。试验结束后应报告实际信息量、效应、CI、事件数与精度。

7 分析集、参与者流程与数据审计

7.1 分析集必须服务于 estimand

分析集 典型定义 主要用途与风险
ITT / Full analysis set 尽可能按随机分配纳入所有参与者 保护随机化;结局缺失仍需合适方法
Modified ITT 增加随机化后条件,例如至少一次给药 若由预后或治疗影响,可能产生选择偏倚
Per protocol 达到预设依从和无重大偏离 不是随机比较;非劣效中常作共同证据
As treated 按实际接受治疗分组 易受换药与依从的混杂,不等同随机效应
Safety set 通常至少接受一次研究治疗,按实际暴露 分母、暴露时间与治疗归类必须明确

ITT 是随机分配策略,不是“用 LOCF 填满数据”。若要估计严格符合方案效应,通常需要处理依从、偏离和随机化后选择,而不只是删掉“违规者”。

7.2 最小审计与流程表

在建模前核对 ID 唯一、分组非缺失、时间顺序、变量范围、主要结局推导、重复记录、分母、暴露、方案偏离、停药与失访原因。下面的表把“未完成治疗”“未测到结局”和“未进入分析”分开;真实报告还需 CONSORT 流程图。

flow_table <- data.frame(
  stage = c("随机分配", "第 12 周结局已测", "第 12 周结局缺失",
            "纳入 NRI 应答分析", "纳入安全性分析", "一年生存结局可用"),
  Control = c(
    sum(trial$arm == "Control"),
    sum(trial$arm == "Control" & !is.na(trial$week12)),
    sum(trial$arm == "Control" & is.na(trial$week12)),
    sum(trial$arm == "Control" & !is.na(trial$responder_nri)),
    sum(trial$arm == "Control" & !is.na(trial$adverse_event)),
    sum(trial$arm == "Control" & !is.na(trial$time_days))
  ),
  Treatment = c(
    sum(trial$arm == "Treatment"),
    sum(trial$arm == "Treatment" & !is.na(trial$week12)),
    sum(trial$arm == "Treatment" & is.na(trial$week12)),
    sum(trial$arm == "Treatment" & !is.na(trial$responder_nri)),
    sum(trial$arm == "Treatment" & !is.na(trial$adverse_event)),
    sum(trial$arm == "Treatment" & !is.na(trial$time_days))
  )
)
knitr::kable(flow_table, caption = "模拟试验的分析分母与随访流程")
模拟试验的分析分母与随访流程
stage Control Treatment
随机分配 178 182
第 12 周结局已测 158 157
第 12 周结局缺失 20 25
纳入 NRI 应答分析 178 182
纳入安全性分析 178 182
一年生存结局可用 178 182

8 结局类型决定效应量与模型

结局 首选描述 常见主模型 至少报告
连续单时点 每组均值/SD、分布 ANCOVA:随访值 ~ 分组 + 基线 + 分层因素 调整均差、95% CI、单位
重复连续 每访视 n、均值/SD MMRM 或合适混合模型/GEE 每时点差、CI、协方差与缺失假设
二元 每组事件数/风险 二项模型、标准化风险;稀疏时精确/惩罚方法 风险、RD、RR;OR 需标明
计数/率 事件数与观察人时 Poisson/负二项,offset 为 log 人时 每组率、IRR、过度离散检查
有序等级 各等级比例 比例优势或预设替代模型 全分布与 common OR,检查假设
时间到事件 KM、生存率、风险人数 log-rank、Cox;非 PH 时 RMST 等 事件数、HR/CI、固定时点绝对风险

“正态/非正态”不是选方法的唯一开关。首先匹配 estimand、随机化单位、结局尺度、随访结构和缺失机制;随后检查模型假设,并用预设替代模型评估稳健性。

9 连续主结局:ANCOVA

9.1 为什么通常调整基线

对有基线测量的连续结局,分析第 12 周值并调整基线,通常比只比较变化值更有效率,也能处理偶然基线差异。预后性基线协变量可提高精度;FDA 2023 协变量调整指南 强调在随机试验中预设使用。随机化分层因素一般也应进入主模型。

fit_ancova <- lm(
  week12 ~ arm + baseline + site + severity,
  data = trial
)
ancova_row <- coef(summary(fit_ancova))["armTreatment", ]
ancova_ci <- confint(fit_ancova, "armTreatment")

change_test <- with(
  subset(trial, !is.na(change)),
  t.test(change[arm == "Treatment"], change[arm == "Control"])
)

continuous_results <- data.frame(
  analysis = c("ANCOVA:第 12 周值 + 基线/分层因素",
               "未调整变化值 Welch 比较"),
  estimate = c(ancova_row["Estimate"], unname(change_test$estimate[1] -
                                                change_test$estimate[2])),
  lower = c(ancova_ci[1], change_test$conf.int[1]),
  upper = c(ancova_ci[2], change_test$conf.int[2]),
  p_value = c(ancova_row["Pr(>|t|)"], change_test$p.value)
)
continuous_results$effect_95_CI <- with(
  continuous_results, fmt_ci(estimate, lower, upper)
)
continuous_results$p_value <- vapply(continuous_results$p_value, format_p,
                                     character(1))
knitr::kable(continuous_results[c("analysis", "effect_95_CI", "p_value")],
             caption = "第 12 周连续结局:Treatment − Control,负值有利")
第 12 周连续结局:Treatment − Control,负值有利
analysis effect_95_CI p_value
Estimate ANCOVA:第 12 周值 + 基线/分层因素 -1.76 (-2.93, -0.59) 0.00342
未调整变化值 Welch 比较 -1.72 (-3.00, -0.44) 0.00863

可测病例 ANCOVA 估计调整均差约为 -1.76 分(95% CI -2.93 至 -0.59)。这是一个有用的模型演示,但由于第 12 周并非人人可测,它不能仅凭“按随机分组建模”就称为完整 ITT 分析;结论还依赖缺失机制,后文将用 MI 与 delta 敏感性分析。

old_par <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 2, 1))
plot(fitted(fit_ancova), residuals(fit_ancova),
     pch = 19, col = grDevices::adjustcolor(trial_palette["treatment"], 0.55),
     xlab = "拟合值", ylab = "残差", main = "方差与函数形式")
abline(h = 0, lty = 2, col = trial_palette["harm"])
qqnorm(residuals(fit_ancova), pch = 19,
       col = grDevices::adjustcolor(trial_palette["navy"], 0.55),
       main = "残差 Q–Q 图")
qqline(residuals(fit_ancova), col = trial_palette["harm"], lwd = 2)
两个诊断图,展示 ANCOVA 残差相对拟合值的散点以及残差正态 Q-Q 图。

ANCOVA 残差诊断:左为拟合值与残差,右为正态 Q–Q 图。

par(old_par)

诊断发现问题时,应考虑稳健标准误、变换、合适的广义模型或预设稳健估计,而不是根据结果不断删异常值。治疗系数是组间调整均值差,不是每位参与者的个体变化。

10 重复连续结局:MMRM 的紧凑应用

重复随访能描述起效时间并利用部分记录。常见 MMRM 把 post-baseline visit 作为分类变量、包含 treatment × visit,并为同一参与者的残差指定灵活协方差;它不需要把缺失值先填成 LOCF。在 MAR 假设下,似然使用已有观测,但 MAR 仍不可由数据“检验为真”。

long_trial <- reshape(
  trial[c("id", "arm", "baseline", "site", "severity",
          "week4", "week8", "week12")],
  varying = c("week4", "week8", "week12"),
  v.names = "score", timevar = "visit_index",
  times = 1:3, direction = "long"
)
long_trial$visit <- factor(
  long_trial$visit_index, levels = 1:3,
  labels = c("Week 4", "Week 8", "Week 12")
)
long_trial <- long_trial[order(long_trial$id, long_trial$visit_index), ]

fit_mmrm <- nlme::gls(
  score ~ baseline + site + severity + visit * arm,
  correlation = nlme::corSymm(form = ~ visit_index | id),
  weights = nlme::varIdent(form = ~ 1 | visit),
  data = long_trial, method = "REML", na.action = na.omit,
  control = nlme::glsControl(opt = "optim")
)

mmrm_beta <- coef(fit_mmrm)
mmrm_vcov <- vcov(fit_mmrm)
contrast_matrix <- matrix(0, 3, length(mmrm_beta),
                          dimnames = list(levels(long_trial$visit),
                                          names(mmrm_beta)))
contrast_matrix[, "armTreatment"] <- 1
contrast_matrix["Week 8", "visitWeek 8:armTreatment"] <- 1
contrast_matrix["Week 12", "visitWeek 12:armTreatment"] <- 1
mmrm_estimate <- drop(contrast_matrix %*% mmrm_beta)
mmrm_se <- sqrt(diag(contrast_matrix %*% mmrm_vcov %*% t(contrast_matrix)))
  # 这里使用大样本正态近似;确证分析应预设小样本自由度方法。
mmrm_table <- data.frame(
  visit = rownames(contrast_matrix),
  estimate = mmrm_estimate,
  lower = mmrm_estimate - qnorm(0.975) * mmrm_se,
  upper = mmrm_estimate + qnorm(0.975) * mmrm_se
)
mmrm_table$effect_95_CI <- with(mmrm_table, fmt_ci(estimate, lower, upper))
knitr::kable(mmrm_table[c("visit", "effect_95_CI")],
             caption = "MMRM 各访视调整均差:Treatment − Control,负值有利")
MMRM 各访视调整均差:Treatment − Control,负值有利
visit effect_95_CI
Week 4 Week 4 -1.57 (-2.57, -0.57)
Week 8 Week 8 -2.15 (-3.22, -1.08)
Week 12 Week 12 -1.75 (-2.92, -0.58)

协方差结构、自由度近似和数值收敛应在 SAP 中规定。若主要 estimand 是第 12 周平均差,其他访视可以解释轨迹,但不能因为某个时间点 p 值更小就改主要结论。更完整的随机效应、GEE 与缺失结构见纵向数据专题。

11 二元结局:绝对风险与相对效应并列

11.1 RD、RR 和 OR 回答不同问题

风险差(RD)直接给出每 100 人多/少多少事件,风险比(RR)给出比例变化,优势比(OR)来自 logistic 模型但在事件常见时可能远离 RR。OR 还具有 non-collapsibility:即使没有混杂,调整与未调整 OR 也可能不同,因此不要把差异自动解释为偏倚被“修复”。

binary_effects <- function(y, arm, conf.level = 0.95) {
  keep <- !is.na(y) & !is.na(arm)
  y <- as.integer(y[keep])
  arm <- droplevels(arm[keep])
  stopifnot(nlevels(arm) == 2L, all(y %in% 0:1))
  z <- qnorm(1 - (1 - conf.level) / 2)
  n0 <- sum(arm == levels(arm)[1]); n1 <- sum(arm == levels(arm)[2])
  e0 <- sum(y[arm == levels(arm)[1]]); e1 <- sum(y[arm == levels(arm)[2]])
  p0 <- e0 / n0; p1 <- e1 / n1
  rd <- p1 - p0
  se_rd <- sqrt(p1 * (1 - p1) / n1 + p0 * (1 - p0) / n0)
  cells <- c(e1, n1 - e1, e0, n0 - e0)
  if (any(cells == 0)) cells <- cells + 0.5
  p1_cc <- cells[1] / sum(cells[1:2])
  p0_cc <- cells[3] / sum(cells[3:4])
  rr <- p1_cc / p0_cc
  se_log_rr <- sqrt(1 / cells[1] - 1 / sum(cells[1:2]) +
                    1 / cells[3] - 1 / sum(cells[3:4]))
  odds_ratio <- (cells[1] * cells[4]) / (cells[2] * cells[3])
  se_log_or <- sqrt(sum(1 / cells))
  data.frame(
    measure = c("风险差", "风险比", "优势比"),
    estimate = c(rd, rr, odds_ratio),
    lower = c(rd - z * se_rd, exp(log(rr) - z * se_log_rr),
              exp(log(odds_ratio) - z * se_log_or)),
    upper = c(rd + z * se_rd, exp(log(rr) + z * se_log_rr),
              exp(log(odds_ratio) + z * se_log_or)),
    control_risk = p0, treatment_risk = p1
  )
}

responder_effects <- binary_effects(trial$responder_nri, trial$arm)
responder_effects$effect_95_CI <- with(
  responder_effects, fmt_ci(estimate, lower, upper, digits = 3)
)
knitr::kable(
  responder_effects[c("measure", "effect_95_CI", "control_risk",
                       "treatment_risk")], digits = 3,
  caption = "缺失按无应答处理后的第 12 周应答效应"
)
缺失按无应答处理后的第 12 周应答效应
measure effect_95_CI control_risk treatment_risk
风险差 0.108 (0.009, 0.207) 0.315 0.423
风险比 1.345 (1.021, 1.771) 0.315 0.423
优势比 1.598 (1.037, 2.461) 0.315 0.423
fit_responder <- glm(
  responder_nri ~ arm + baseline + site + severity,
  family = binomial(), data = trial
)
adjusted_or <- exp(c(
  estimate = coef(fit_responder)["armTreatment"],
  confint.default(fit_responder, "armTreatment")
))
data.frame(
  measure = "调整 OR",
  estimate = adjusted_or[1], lower = adjusted_or[2], upper = adjusted_or[3]
) |>
  knitr::kable(digits = 3, caption = "预设协变量调整的 logistic 模型")
预设协变量调整的 logistic 模型
measure estimate lower upper
estimate.armTreatment 调整 OR 1.6 1.04 2.48

缺失按无应答(NRI)把缺失与失败合成,是一个 composite 策略示例,并非普遍“保守”:若两组缺失率或缺失原因不同,方向可能复杂。正式分析应同时给每组原始分子/分母、绝对效应和敏感性分析。稀疏事件或完全分离时,应考虑精确或惩罚方法,不能依赖普通 Wald CI。

12 计数与发生率:把观察时间放进 offset

fit_poisson <- glm(
  episode_count ~ arm + baseline + site + severity +
    offset(log(exposure_weeks)),
  family = poisson(), data = trial
)
dispersion <- sum(residuals(fit_poisson, type = "pearson")^2) /
  df.residual(fit_poisson)
fit_quasipoisson <- update(fit_poisson, family = quasipoisson())
count_row <- coef(summary(fit_quasipoisson))["armTreatment", ]
irr <- exp(c(
  estimate = count_row["Estimate"],
  lower = count_row["Estimate"] - qnorm(0.975) * count_row["Std. Error"],
  upper = count_row["Estimate"] + qnorm(0.975) * count_row["Std. Error"]
))
rate_table <- aggregate(
  cbind(events = episode_count, person_weeks = exposure_weeks) ~ arm,
  data = trial, FUN = sum
)
rate_table$rate_per_100_person_weeks <-
  100 * rate_table$events / rate_table$person_weeks
knitr::kable(rate_table, digits = 2,
             caption = "两组急性发作数、观察人周与粗发生率")
两组急性发作数、观察人周与粗发生率
arm events person_weeks rate_per_100_person_weeks
Control 271 1784 15.2
Treatment 207 1806 11.5
knitr::kable(data.frame(
  IRR = irr[1], lower = irr[2], upper = irr[3],
  Pearson_dispersion = dispersion
), digits = 3, caption = "准 Poisson 模型的 Treatment/Control IRR")
准 Poisson 模型的 Treatment/Control IRR
IRR lower upper Pearson_dispersion
estimate.Estimate 0.746 0.587 0.948 1.75

offset 的系数固定为 1,把不同观察时间转换为率。Pearson 离散度明显大于 1 时,普通 Poisson 标准误偏小;准 Poisson 可调整精度但不改变点估计,负二项模型还改变均值–方差关系。复发事件若存在个体内相关、终末事件或事件依赖,还应考虑 GEE、脆弱性或专门复发事件模型。

13 时间到事件结局

13.1 KM、log-rank 与 Cox 各有角色

Kaplan–Meier 估计随时间的生存概率,log-rank 检验总体曲线差异,Cox 模型估计条件风险比。HR 不是风险比,也不是“复发时间延长 34%”。删失必须在给定协变量与模型假设下具有可辩护的独立性。

survival_outcome <- with(trial, survival::Surv(time_days, status))
fit_km <- survival::survfit(survival_outcome ~ arm, data = trial)
fit_logrank <- survival::survdiff(survival_outcome ~ arm, data = trial)
logrank_p <- pchisq(fit_logrank$chisq,
                    df = length(fit_logrank$n) - 1, lower.tail = FALSE)
fit_cox <- survival::coxph(
  survival_outcome ~ arm + baseline + site + severity,
  data = trial
)
cox_row <- summary(fit_cox)$coefficients["armTreatment", ]
cox_ci <- exp(c(coef(fit_cox)["armTreatment"],
                confint(fit_cox, "armTreatment")))

km_summary <- summary(fit_km, times = c(90, 180, 365), extend = TRUE)
km_table <- data.frame(
  arm = sub("^arm=", "", km_summary$strata),
  day = km_summary$time,
  survival = km_summary$surv,
  lower = km_summary$lower,
  upper = km_summary$upper
)
knitr::kable(km_table, digits = 3,
             caption = "第 90、180 与 365 天 KM 无复发生存概率")
第 90、180 与 365 天 KM 无复发生存概率
arm day survival lower upper
Control 90 0.857 0.806 0.910
Control 180 0.696 0.630 0.769
Control 365 0.485 0.413 0.568
Treatment 90 0.892 0.848 0.939
Treatment 180 0.782 0.721 0.847
Treatment 365 0.615 0.543 0.696
knitr::kable(data.frame(
  analysis = c("Log-rank", "调整 Cox HR"),
  estimate = c(NA, cox_ci[1]),
  lower = c(NA, cox_ci[2]), upper = c(NA, cox_ci[3]),
  p_value = c(logrank_p, cox_row["Pr(>|z|)"])
), digits = 3, caption = "一年内复发/进展的组间比较")
一年内复发/进展的组间比较
analysis estimate lower upper p_value
Log-rank NA NA NA 0.014
armTreatment 调整 Cox HR 0.658 0.474 0.913 0.012
plot(fit_km, col = c(trial_palette["control"], trial_palette["treatment"]),
     lwd = 2, mark.time = TRUE, conf.int = FALSE,
     xlab = "随机化后天数", ylab = "无复发生存概率", ylim = c(0.35, 1))
legend("bottomleft", legend = c("Control", "Treatment"),
       col = c(trial_palette["control"], trial_palette["treatment"]),
       lwd = 2, bty = "n")
Control 与 Treatment 两条 Kaplan-Meier 无复发生存曲线,Treatment 曲线总体较高。

模拟试验的一年 Kaplan–Meier 无复发生存曲线;短竖线表示删失。

ph_check <- survival::cox.zph(fit_cox)
knitr::kable(as.data.frame(ph_check$table), digits = 3,
             caption = "Schoenfeld 残差比例风险检验;低把握度下不能只凭 p 值")
Schoenfeld 残差比例风险检验;低把握度下不能只凭 p 值
chisq df p
arm 0.000 1 0.998
baseline 0.252 1 0.616
site 4.198 3 0.241
severity 0.688 1 0.407
GLOBAL 4.885 6 0.559

比例风险不合理时,预设固定时间生存差、限制平均生存时间(RMST)、分段效应或时间变化系数往往更直接。竞争事件存在时,需区分原因别风险与累积发生风险;把竞争事件当普通删失不一定回答目标问题。详见生存分析专题。

14 缺失数据与敏感性分析

14.1 先预防,再建模,最后挑战假设

应优先减少缺失:即使停止治疗,也尽量继续收集与 estimand 有关的结局和 intercurrent-event 时间、原因。完整病例分析无偏需要强条件;LOCF 通常扭曲轨迹与方差;“MAR”不是一种填补算法,而是给定已观测信息后缺失不再依赖未观测值的假设。

EMA 缺失数据指南 强调预防、主分析假设和敏感性。MI 模型应包含结局预测因子、缺失预测因子、分组、分层变量及分析模型的重要结构,并按 Rubin 规则合并;一次均值填补会虚假提高精度。

14.2 MAR 多重插补与 delta pattern-mixture

下面先在每组内以已观测第 12 周结果拟合插补模型;再把治疗组插补值统一增加 δ\delta 分,表示“在相同已观测信息下,治疗组缺失参与者比 MAR 预测更差”。每个 δ\delta 使用相同随机数(common random numbers),让变化主要来自假设而非 Monte Carlo 噪声。

mi_ancova <- function(data, delta = 0, m = 100L, seed = 41000L) {
  set.seed(seed)
  estimates <- variances <- numeric(m)
  for (j in seq_len(m)) {
    completed <- data
    for (group in levels(completed$arm)) {
      observed <- completed$arm == group & !is.na(completed$week12)
      missing <- completed$arm == group & is.na(completed$week12)
      if (!any(missing)) next
      imputation_fit <- lm(
        week12 ~ baseline + site + severity + sex,
        data = completed, subset = observed
      )
      beta_hat <- coef(imputation_fit)
      beta_vcov <- vcov(imputation_fit)
      beta_draw <- drop(beta_hat +
                          t(chol(beta_vcov)) %*% rnorm(length(beta_hat)))
      x_missing <- model.matrix(
        ~ baseline + site + severity + sex,
        data = completed[missing, , drop = FALSE]
      )
      imputed <- drop(x_missing %*% beta_draw) +
        rnorm(sum(missing), 0, summary(imputation_fit)$sigma)
      if (group == "Treatment") imputed <- imputed + delta
      completed$week12[missing] <- imputed
    }
    analysis_fit <- lm(
      week12 ~ arm + baseline + site + severity,
      data = completed
    )
    estimates[j] <- coef(analysis_fit)["armTreatment"]
    variances[j] <- vcov(analysis_fit)["armTreatment", "armTreatment"]
  }
  q_bar <- mean(estimates)
  u_bar <- mean(variances)
  between <- var(estimates)
  total <- u_bar + (1 + 1 / m) * between
  df <- if (between < .Machine$double.eps) Inf else
    (m - 1) * (1 + u_bar / ((1 + 1 / m) * between))^2
  critical <- qt(0.975, df)
  data.frame(
    delta = delta, estimate = q_bar, se = sqrt(total), df = df,
    lower = q_bar - critical * sqrt(total),
    upper = q_bar + critical * sqrt(total),
    p_value = 2 * pt(-abs(q_bar / sqrt(total)), df)
  )
}

delta_grid <- seq(0, 8, by = 0.5)
tipping <- do.call(rbind, lapply(delta_grid, function(delta_i) {
  mi_ancova(trial, delta = delta_i, m = 100L, seed = 41000L)
}))

selected_tipping <- tipping[tipping$delta %in% c(0, 2, 4, 5, 5.5, 6, 8), ]
selected_tipping$effect_95_CI <- with(
  selected_tipping, fmt_ci(estimate, lower, upper, digits = 3)
)
selected_tipping$p_value <- vapply(selected_tipping$p_value, format_p,
                                   character(1))
knitr::kable(selected_tipping[c("delta", "effect_95_CI", "p_value")],
             caption = "治疗组缺失值 delta 调差后的 MI 敏感性分析")
治疗组缺失值 delta 调差后的 MI 敏感性分析
delta effect_95_CI p_value
1 0.0 -1.877 (-3.047, -0.708) 0.00166
5 2.0 -1.605 (-2.777, -0.434) 0.00726
9 4.0 -1.333 (-2.516, -0.151) 0.0271
11 5.0 -1.197 (-2.389, -0.006) 0.0489
12 5.5 -1.129 (-2.326, 0.067) 0.0643
13 6.0 -1.061 (-2.263, 0.141) 0.0835
17 8.0 -0.789 (-2.019, 0.440) 0.208
plot(tipping$delta, tipping$estimate, type = "b", pch = 19,
     col = trial_palette["treatment"], lwd = 2,
     ylim = range(tipping$lower, tipping$upper, 0),
     xlab = expression(paste("治疗组缺失值调差 ", delta, "(分)")),
     ylab = "Treatment − Control 调整均差")
segments(tipping$delta, tipping$lower, tipping$delta, tipping$upper,
         col = grDevices::adjustcolor(trial_palette["treatment"], 0.55))
abline(h = 0, lty = 2, col = trial_palette["harm"])
治疗效应及95%置信区间随治疗组缺失结局delta调差而向零移动,在delta约5.5时置信区间跨过零。

Delta-adjusted pattern-mixture 敏感性曲线;横线为无效值 0。

MAR(δ=0\delta=0)结果约为 -1.88 分;直到把每个治疗组缺失结局平均调差约 5.5 分,双侧 95% CI 才首次包含 0。tipping point 的意义取决于 5.5 分是否临床可信,不是简单宣布“稳健”。真实试验还应让 delta 与停药原因、访视时间和专家知识对应,并检查插补模型与主分析相容性。

15 多重性与期中分析

15.1 多个终点不是三个独立故事

多个主要终点、多个剂量、多个时间点和反复查看数据都会增加错误阳性。优先顺序/层级检验、gatekeeping、Holm 或 Bonferroni 控制家族错误率(FWER);FDR 更常用于探索性大量假设。调整方案应在 SAP 中定义假设家族,而不是把所有 p 值机械塞进同一列表。

p_continuous <- coef(summary(fit_ancova))["armTreatment", "Pr(>|t|)"]
p_responder <- coef(summary(fit_responder))["armTreatment", "Pr(>|z|)"]
p_survival <- summary(fit_cox)$coefficients["armTreatment", "Pr(>|z|)"]
multiplicity <- data.frame(
  endpoint = c("第 12 周连续评分", "第 12 周应答", "一年复发/进展"),
  raw_p = c(p_continuous, p_responder, p_survival)
)
multiplicity$holm_p <- p.adjust(multiplicity$raw_p, method = "holm")
multiplicity$reject_at_0.05 <- multiplicity$holm_p < 0.05
knitr::kable(multiplicity, digits = 4,
             caption = "三个假设被预先定义为同一家族时的 Holm 调整示例")
三个假设被预先定义为同一家族时的 Holm 调整示例
endpoint raw_p holm_p reject_at_0.05
第 12 周连续评分 0.0034 0.0103 TRUE
第 12 周应答 0.0327 0.0327 TRUE
一年复发/进展 0.0124 0.0248 TRUE

15.2 期中查看必须花费 alpha

独立 DMC 按章程审阅疗效、无效和安全性;发起人运营团队通常不应看到非盲趋势。信息比例应尽可能基于事件数或估计信息,而不只是日历时间。下面只把累计 alpha 转换成“等效双侧 z 阈值”,用于理解早期门槛为何严格;它不是考虑各次检验相关性的精确序贯边界。

alpha <- 0.05
information_fraction <- c(0.33, 0.67, 1.00)
z_alpha <- qnorm(1 - alpha / 2)
spend_obf <- 2 - 2 * pnorm(z_alpha / sqrt(information_fraction))
spend_pocock <- alpha * log(1 + (exp(1) - 1) * information_fraction)
spending <- rbind(
  data.frame(method = "O'Brien–Fleming-like", look = 1:3,
             information_fraction = information_fraction,
             cumulative_alpha = spend_obf),
  data.frame(method = "Pocock-like", look = 1:3,
             information_fraction = information_fraction,
             cumulative_alpha = spend_pocock)
)
spending$incremental_alpha <- ave(
  spending$cumulative_alpha, spending$method,
  FUN = function(x) c(x[1], diff(x))
)
spending$cumulative_z_equivalent <-
  qnorm(1 - spending$cumulative_alpha / 2)
knitr::kable(spending, digits = 4,
             caption = "教学用累计 alpha spending;不是精确相关序贯边界")
教学用累计 alpha spending;不是精确相关序贯边界
method look information_fraction cumulative_alpha incremental_alpha cumulative_z_equivalent
O’Brien–Fleming-like 1 0.33 0.0006 0.0006 3.41
O’Brien–Fleming-like 2 0.67 0.0166 0.0160 2.39
O’Brien–Fleming-like 3 1.00 0.0500 0.0334 1.96
Pocock-like 1 0.33 0.0225 0.0225 2.28
Pocock-like 2 0.67 0.0383 0.0158 2.07
Pocock-like 3 1.00 0.0500 0.0117 1.96

确证设计应使用经过验证的软件计算相关边界、条件把握度和 operating characteristics,并预设非约束/约束无效边界、样本量重估和过度运行数据的处理。安全原因可随时行动,但疗效提前停止仍需校正估计与报告信息时间。

16 亚组分析:直接检验交互

fit_interaction <- lm(
  week12 ~ arm * sex + baseline + site + severity,
  data = trial
)
interaction_beta <- coef(fit_interaction)
interaction_vcov <- vcov(fit_interaction)
interaction_df <- df.residual(fit_interaction)
main_term <- "armTreatment"
interaction_term <- "armTreatment:sexMale"

female_est <- interaction_beta[main_term]
female_se <- sqrt(interaction_vcov[main_term, main_term])
male_est <- interaction_beta[main_term] + interaction_beta[interaction_term]
male_se <- sqrt(interaction_vcov[main_term, main_term] +
                  interaction_vcov[interaction_term, interaction_term] +
                  2 * interaction_vcov[main_term, interaction_term])
critical <- qt(0.975, interaction_df)
subgroup <- data.frame(
  subgroup = c("Female", "Male"),
  estimate = c(female_est, male_est),
  lower = c(female_est - critical * female_se,
            male_est - critical * male_se),
  upper = c(female_est + critical * female_se,
            male_est + critical * male_se)
)
interaction_p <- coef(summary(fit_interaction))[
  interaction_term, "Pr(>|t|)"
]
subgroup$effect_95_CI <- with(subgroup, fmt_ci(estimate, lower, upper))
knitr::kable(subgroup[c("subgroup", "effect_95_CI")],
             caption = paste0("探索性性别亚组;交互 p = ", format_p(interaction_p)))
探索性性别亚组;交互 p = 0.432
subgroup effect_95_CI
Female -1.37 (-2.91, 0.17)
Male -2.33 (-4.16, -0.50)
plot(subgroup$estimate, seq_len(nrow(subgroup)),
     xlim = range(subgroup$lower, subgroup$upper, 0),
     ylim = c(0.5, nrow(subgroup) + 0.5), yaxt = "n",
     pch = 19, col = trial_palette["treatment"],
     xlab = "调整均差(Treatment − Control)", ylab = "")
segments(subgroup$lower, seq_len(nrow(subgroup)),
         subgroup$upper, seq_len(nrow(subgroup)),
         col = trial_palette["treatment"], lwd = 2)
axis(2, at = seq_len(nrow(subgroup)), labels = subgroup$subgroup, las = 1)
abline(v = 0, lty = 2, col = trial_palette["harm"])
Female 和 Male 两个亚组的Treatment减Control均差及95%置信区间;一组区间跨零,另一组未跨零,但交互检验不显著。

探索性性别亚组的调整均差;虚线为无效值 0。

Female 的 CI 包含 0,Male 的 CI 不包含 0,但 treatment × sex 交互 p 约为 0.43:“一组显著、另一组不显著”不是组间差异显著。亚组应少量预设、直接检验交互、优先保留连续修饰变量,并结合方向、可信度、多重性和生物机制解释。EMA 亚组指南 可作为进一步阅读。

17 安全性分析与获益–风险

ae_effects <- binary_effects(trial$adverse_event, trial$arm)
ae_rr <- ae_effects[ae_effects$measure == "风险比", ]
ae_fisher <- fisher.test(table(trial$arm, trial$adverse_event))
ae_summary <- data.frame(
  arm = levels(trial$arm),
  participants = as.integer(table(trial$arm)),
  with_event = as.integer(tapply(trial$adverse_event, trial$arm, sum)),
  risk = as.numeric(tapply(trial$adverse_event, trial$arm, mean))
)
knitr::kable(ae_summary, digits = 3,
             caption = "至少一次模拟不良事件的参与者")
至少一次模拟不良事件的参与者
arm participants with_event risk
Control 178 20 0.112
Treatment 182 38 0.209
knitr::kable(data.frame(
  measure = "Treatment/Control RR",
  estimate = ae_rr$estimate, lower = ae_rr$lower, upper = ae_rr$upper,
  Fisher_p = ae_fisher$p.value
), digits = 3, caption = "不良事件相对风险与精确检验")
不良事件相对风险与精确检验
measure estimate lower upper Fisher_p
Treatment/Control RR 1.86 1.13 3.06 0.015

安全性分析应说明 safety set、实际暴露、风险窗口、事件严重程度、严重不良事件(SAE)、停药、死亡和因果关联评价。常见事件可报告风险差/比,暴露时间不等时报告率;罕见严重伤害通常以 CI、累积病例、叙述审阅和跨来源证据为主。对几十种 AE 各做未校正 p 值并筛“显著项”,既缺乏把握度也误导获益–风险判断。

18 从 SAP 到 CONSORT 2025 报告

18.1 一个可执行的 SAP 最小清单

  • 版本、日期、揭盲状态和偏离处理;
  • 每个主要/关键次要 estimand 的五个属性;
  • 数据快照、分析集、变量推导、访视窗与基线定义;
  • 随机化分层因素、协变量、模型、对比、效应方向与 CI;
  • ICE、缺失数据主假设、敏感性和补充分析;
  • 多重性、期中分析、亚组、安全性、模型失败与软件版本;
  • 预设表图和可追溯的质量控制/独立复核。

18.2 报告效应,不只报告 p 值

一个合格结果句可以写成:

在有第 12 周测量的模拟参与者中,按随机分组并调整基线、中心和基线严重程度后,Treatment 相对 Control 的症状评分平均低 1.76 分(95% CI 0.59 至 2.93 分更低;双侧 p=0.003)。这一 available-case 估计依赖缺失机制;MAR MI 和预设 MNAR delta 分析用于评估稳健性。

结果段还应给两组分母与描述值、实际随访、事件数、绝对效应、模型和方向。统计显著不等于临床重要;CI 同时表达可兼容的获益和伤害范围。

CONSORT 2025 已取代 CONSORT 2010,核心是 30 项清单和参与者流程图,并强化注册、protocol/SAP 可访问、数据与代码共享、患者/公众参与、伤害、实际干预实施、分析集和缺失数据。报告指南提高透明度,但不能修复差的设计。

19 常见错误与修正

常见错误 为什么错 更好的做法
基线逐项做 p 值,显著才调整 随机差异不应驱动模型 按预后价值、分层和 SAP 预设
只做组内前后检验 组内显著性不等于组间差异 直接估计组间效应与 CI
ITT = LOCF 分析策略与填补法被混为一谈 先定义 estimand,再选择缺失方法
p>0.05 就“等效” 缺少精度不证明相似 用预设界值和适当单/双侧 CI
一亚组显著、另一亚组不显著 没有检验效应差异 直接检验 treatment × subgroup
HR 当作风险比 瞬时条件效应不同于累计风险 同报 KM 固定时点风险;检查 PH
每周查看 p 值,显著即停 I 类错误和估计偏倚增加 预设边界、DMC 与 alpha spending
只报 OR 和 p 值 临床绝对意义不清 同报每组风险、RD/RR 与 CI
排除停药者做“干净 ITT” 破坏随机化并选择人群 FAS 主分析;PP 作预设补充/因果方法
安全性逐项找显著 低把握度与多重性严重 描述分母、暴露、严重性、CI 与病例

知识检查与练习

  1. 停药后仍收集第 12 周评分,属于哪种 estimand 策略?

若分析无论停药都使用该评分,通常对应 treatment-policy 策略。停药是 ICE,评分是否缺失是另一个问题。

  1. 为什么分配隐藏和盲法不能互换?

分配隐藏防止随机化前的选择性入组;盲法降低随机化后行为、共同干预和结局评价受知晓分组影响。两者发生时间和针对偏倚不同。

  1. 非劣效试验 p>0.05,能否得出非劣效?

不能。必须检查方向一致的单侧置信上界是否低于预设且有临床/历史依据的界值,并评价 assay sensitivity、FAS/PP 与偏离。

  1. 为什么 ANCOVA 结果不能自动称为完整 ITT?

模型按随机组比较,但只纳入有第 12 周值者。缺失参与者的结局仍需与 estimand 一致的方法和敏感性分析。

  1. OR=1.60 是否表示应答概率增加 60%?

不是。OR 比较 odds;本例应同时查看原始风险、RD 与 RR。事件常见时 OR 与 RR 可明显不同。

  1. Female 不显著、Male 显著,是否证明性别修饰疗效?

不证明。本例直接交互检验 p≈0.43;应解释交互估计与 CI,并考虑探索性和多重性。

  1. Delta tipping point 是“通过/失败”阈值吗?

不是。它说明结论需要多大未观测偏离才改变;关键是该 delta 与临床知识、缺失原因及数据范围相比是否可信。

  1. 为什么期中分析不能把三次 p 值都与 0.05 比?

反复机会增加家族 I 类错误。需要预设相关序贯边界或 alpha spending,并由独立 DMC 按章程执行。

  1. 设计一个自己的主要 estimand。

依次写出 population、treatment conditions、variable、每种 ICE 的策略和 population-level summary,再说明相应主 estimator 与至少一个关键敏感性分析。

快速参考与最终检查清单

分析前

  • 临床问题、效应方向、主要时间点和 MCID 已明确;
  • estimand 五要素与 ICE 策略逐项写清;
  • 随机化、隐藏、盲法、紧急揭盲和审计责任可执行;
  • 样本量考虑现实参数、失访、设计效应、多重性和期中计划;
  • protocol、注册和 SAP 在非盲结果可见前完成并版本化。

数据与模型

  • ID、分母、访视、暴露、事件、偏离和缺失原因已核对;
  • 模型匹配结局尺度、随机化单位、分层与重复结构;
  • 每个结果给效应量、95% CI、单位、方向和每组描述值;
  • 缺失主假设与至少一个针对不可验证假设的敏感性分析对应;
  • 多重性、期中、亚组与安全性未按结果临时改规则;
  • 模型诊断、收敛、替代模型与软件版本留有审计记录。

报告与解释

  • 按 CONSORT 2025 给参与者流程、分配、实际干预、分析人数和伤害;
  • 报告 protocol/SAP/注册号、修订、代码与可共享数据路径;
  • 区分统计显著、临床重要、不精确和无证据;
  • 不把 OR 当 RR、HR 当累计风险比、高 p 值当等效;
  • 结论局限在预设人群、干预、对照、结局、时间和实施环境。
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    xfun_0.60      
##  [6] lattice_0.22-9  splines_4.6.1   cachem_1.1.0    knitr_1.51      htmltools_0.5.9
## [11] rmarkdown_2.31  stats4_4.6.1    lifecycle_1.0.5 cli_3.6.6       grid_4.6.1     
## [16] sass_0.4.10     jquerylib_0.1.4 compiler_4.6.1  tools_4.6.1     nlme_3.1-169   
## [21] evaluate_1.0.5  bslib_0.12.0    survival_3.8-6  yaml_2.3.12     rlang_1.3.0    
## [26] jsonlite_2.0.0  MASS_7.3-65