手里攥着几十万条观察性记录被追问这两种起始策略到底哪个更好的时候真正的难点从来不是跑不出一个显著的数字而是这个数字到底能不能当因果效应来读。目标模拟试验Target Trial EmulationTTE就是冲着这件事来的框架它要求你先像设计一项随机对照试验那样把入组、分配、随访起点、结局、因果对比全部写死然后再回到观察性数据里去拼时间轴、造克隆、加权重把人为引入的偏差一层层掰回来。这套方法在流行病学和药物效果比较里已经用了不少年但真正做到照着流程走一遍就能出可信结论的人并不多大部分翻车都发生在几个不起眼的时间点和权重细节上。这篇东西写给两类人一类是手里有队列数据、理赔数据或者电子病历被业务方或导师要求给个因果结论的分析岗另一类是想把观察性研究做得更硬一点、不想每次都被审稿人问你怎么排除混杂的研究者。我会按实际动手的顺序讲先搞清楚这框架到底在补什么缺口再把研究问题写成一份可执行的方案然后是克隆、删失、加权这套核心机器怎么转估计量怎么选、方差怎么算最后是我自己踩过的几个坑。全程会给出可以直接抄的代码片段和判断阈值但不会替你决定研究问题——那部分只能你自己想清楚。1. 观察性数据想冒充随机试验缺的到底是什么1.1 从组间差异显著到因果效应可解释之间的三道坎直接比较两组人的结局得到的那个数字在同一批数据里几乎必然是有偏的而且偏的方向你还未必知道。第一道坎是混杂谁用 A 方案、谁用 B 方案从来不是抛硬币决定的。临床实践中病情更重的往往先上强方案经济条件更好、依从性更高的人更容易拿到新药这些因素同时影响用哪个方案和结局好坏把它们混在一起算出来的差异本质上是人群构成的差异不是方案本身的差异。第二道坎是时间零点错位也就是不死时间偏倚。举个最常见的错误做法把确诊日定为随访起点把是否用过某方案作为分组依据。问题在于从确诊到用上药之间有一段时间这段时间里死亡或者病情急转直下的人永远进不了用药组。于是用药组天然就活得久一点这个优势是筛选出来的跟药本身没关系。我见过的返修意见里这条被点名的频率最高。第三道坎是依从性与方案偏离。真实世界里的人会换药、会停药、会因为副作用改成隔天吃。如果你按实际吃了什么分组就等于在比较一群能坚持下来的人和一群坚持不下来的人而能坚持本身就是预后更好的代理指标。反过来按起始方案分组又没法回答如果所有人都严格遵守会怎样。这两类问题在随机试验里靠 ITT 和 per-protocol 两套分析解决观察性数据里没有随机化就得靠 TTE 这一套机制去模拟。1.2 TTE 真正改变的是工作顺序不是统计方法很多人以为 TTE 是个新估计量其实它不是。它用的是加权、标准化、pooled logistic 这些早就有的工具真正被改变的是做事的顺序。传统分析路径是先看看数据里有什么跑个回归发现系数不对劲再回头补几个敏感性分析最后在讨论部分写一段可能存在残余混杂。TTE 把这条路反过来先写方案方案里每个组件都必须能在数据里对应到具体的一列或者一个判定规则写不出来的组件说明你的数据根本支撑不了这个问题。我习惯用一个很粗暴的自检如果你没法用一句话说清楚时间零点是哪一天这个研究就别做了。因为时间零点定义不清楚后面的随访、删失、加权全部失去锚点做出来的东西没法解释。这个判断标准看起来武断但在实际项目里救过我很多次——有好几次需求方想要的效果数据里的时间戳精度只到月份那这个研究从根上就不成立。顺序反过来的好处是很实在的你在看结局数据之前就把所有决策固定下来后面不管跑出什么结果都没法微调。这跟预注册是一个道理只不过 TTE 把它变成了流程的一部分而不只是一份伦理要求。1.3 这套框架能接什么活接不了什么活TTE 不是万能的把边界说清楚比吹它的能力更有用。它成立依赖三个前提任何一个破了结论就得打折扣。前提具体含义破了会怎样可交换性在给定基线协变量后两组人的潜在结局分布可比残余混杂效应估计偏离真值正值性每个策略在每个协变量层里都有人执行权重爆炸估计量方差失控一致性策略定义明确观察到的结局确实对应所定义的策略策略版本混杂效应含义模糊正值性这条最容易被忽视。假如 B 方案在某个年龄层里几乎没人用你在那一层就没有可比对象加权会把少数几个人的权重放大到几十甚至上百最终估计值被这几个人绑架。一致性也一样如果A 方案在数据里既包含标准剂量也包含减量那你估出来的就不是起始 A 方案的效应而是这些版本的加权平均解释起来非常别扭。还有一类问题 TTE 完全帮不上忙测量误差。如果结局判定本身就有系统性漏诊或者关键混杂变量根本没记录那再精巧的设计也只是在错误的地基上盖楼。我一般会在项目启动会上就把这句话摆出来省得后面扯皮。2. 把研究问题钉成一份能执行的试验方案2.1 八个组件逐条落地一条都不能空着TTE 的方案文档通常拆成八个组件。这八个不是学术摆设每一个都直接对应后面的一段代码或者一个数据操作。我把它整理成一张对照表左边是随机试验里的样子右边是观察性数据里的落地方式中间那列是你在写方案时必须回答的问题。组件随机试验里的形态观察性数据落地要点入组标准筛选合格受试者所有筛选变量必须在时间零点当天或之前测量完成不能用未来信息治疗策略随机分配 A 或 B明确定义起始 A与起始 B并说明后续换药是否允许分配程序随机化观察到的实际起始方案需声明无未测量混杂假设时间零点随机化当天入组合格、起始用药、随访开始三者必须是同一天结局预定义终点预定义判定规则、判定窗口与数据来源随访期从随机化到事件或截止从时间零点到结局、删失或行政截尾因果对比ITT 或 per-protocol与所选估计目标一一对应不能混着说分析计划预先写定加权方式、截尾位置、方差方法、敏感性分析全部预先固定实际写的时候最容易空着的是治疗策略和时间零点这两栏。策略定义得好不好直接决定了后面克隆那一步能不能做。我的经验是策略要写成在时间零点开始使用 X 类药物且在随访期内不切换到其他类别这种可判定的句子而不是使用 X 治疗这种含糊表述。前者在数据里能落到一行具体的处方记录后者你没法写代码去判。参数上还有几个要提前拍板的宽限期给多长常见是 30 天、60 天或 90 天取决于临床上开始用药的实际延迟分布、随访上限定多长受数据可得时间限制但也要考虑临床意义、结局判定窗口多宽比如诊断编码出现一次就算还是要求两次或住院记录佐证。这些数字不是拍脑袋最好画一下用药起始时间的分布图看看 90% 的人落在多少天以内宽限期就设在那个位置附近。2.2 时间零点对齐整个框架最脆的一环时间零点必须同时满足三件事个体符合入组标准、被分配到某个策略、随访开始计时。三者在观察性数据里天然是错开的把它们强行对齐是 TTE 最费功夫的一步。实际操作中我一般按这个顺序处理先确定起始用药事件第一条处方记录或第一次给药记录把它作为候选时间零点然后回溯检查这个时点上所有入组标准是否成立如果某个人的入组标准依赖于起始用药之后才能获得的信息比如用药后三个月内做过某检查这个人必须被剔除因为这是对撞变量用它筛人会引入选择偏倚。这一点反直觉很多人在写纳排标准时会不自觉地写进去。# 伪代码把起始用药日作为候选零点并回溯校验入组标准 library(data.table) base - dt[order(id, rx_date)][, .SD[1], by id] # 每个个体的首次用药 base - base[, t0 : rx_date] # 候选时间零点 # 入组标准必须全部在 t0 当天或之前成立 base - base[ age_at_t0 18 prior_days_enrolled 180 # 有足够的基线观察窗 is.na(pregnancy_flag_at_t0) !(dx_code %in% exclusion_codes_before_t0) ] # 剔除在 t0 之后才产生信息的筛选条件所导致的个体 base - base[!(id %in% ids_with_post_t0_criteria)]这段代码看着简单但prior_days_enrolled 180这一行是有讲究的。基线观察窗太短你无法确认这个人真的是新使用者可能之前在其他地方已经用过同类药物窗太长样本量掉得厉害。常见折中是 180 天或 365 天具体选哪个看你数据覆盖的连续性和目标人群的用药周期。我一般会两个都跑一遍看结论是否稳健如果差很多说明基线定义对结果影响大需要在正文里交代。新使用者设计配合活性对照是这套框架的标配组合。活性对照的意思是两组都用某种真实存在的治疗而不是一组用药一组什么都不做。后者的问题在于什么都不做往往对应着一群更健康或者更没机会看病的人混杂结构完全不可控。这个选择在方案阶段就要定下来不能等到跑完发现结果不好看再换。2.3 方案固化在看结局数据之前把门关死写方案这件事形式上是一份文档实质上是一次自我约束。我的做法是建一个版本化的方案文件包含所有决策点、决策理由、以及日期戳。所有参与的人在上面签字哪怕是口头确认也行然后才开始碰结局数据。具体要固定下来的东西包括纳排标准的确切定义和编码列表、时间零点的判定规则、宽限期长度、随访上限、结局定义与判定窗口、删失规则什么情况算删失、什么情况算竞争事件、权重模型里放哪些协变量以及为什么放这些、权重截尾的分位数、方差估计方法、预先计划的敏感性分析清单。这份清单看起来很长但真正写下来也就两三页而且写的过程中经常能发现自己思路里的漏洞。有个细节值得单独说协变量清单要在看结局之前定不能靠单变量筛选或者逐步回归来挑。后者的逻辑是和结局相关的才放进模型这在预测建模里合理在因果推断里是错的因为混杂变量的选择依据应该是它是否同时影响处理和结局跟它在样本里是否显著没关系。一个弱效应但确实是混杂的变量被逐步回归踢掉了偏倚就留下了。我一般按领域知识列一份清单宁可多放几个然后用有效样本量去检查权重会不会被撑爆。3. 克隆、删失、加权让观察性数据长得像随机化3.1 克隆这一步是在给每个人补上另一条路克隆法的核心思想很直白既然随机试验里每个人都有可能被分到 A 或 B那观察性数据里我也给每个人造出这两个版本。在时间零点把每个个体复制成 N 份N 等于策略数每份贴上不同的策略标签然后各自按自己的策略往下走。# 每个个体在 t0 克隆成两份分别赋策略 A / B clone_A - base[, .(id, clone_id paste0(id, _A), strategy A, t0, age_at_t0, sex, comorbidity_score, index_year)] clone_B - base[, .(id, clone_id paste0(id, _B), strategy B, t0, age_at_t0, sex, comorbidity_score, index_year)] clones - rbind(clone_A, clone_B)关键在于接下来这一步每个克隆只有在其实际数据与该策略一致的时候才继续被随访。也就是说一个实际用了 A 的人他的 A 克隆会一直跟到结局而他的 B 克隆在时间零点就被人为删失了因为他没有起始 B。反过来一个实际用了 B 的人B 克隆继续A 克隆立刻删失。这样一来每个策略组里既有实际执行了该策略的人的完整轨迹也有实际没执行但理论上可以执行的人被删失掉。前者提供信息后者被删失这件事本身携带了信息——它告诉我们这个人按协变量分布本来有机会执行这个策略。加权就是用来把这部分信息折算回去的。新手最容易搞混的一点是克隆不是抽样不会增加信息量。它只是把反事实这件事在数据结构上显式表达出来真正干活的是后面的加权。所以克隆完之后样本行数翻倍但有效样本量并不会翻倍这个直觉要建立起来不然看到权重诊断时会一头雾水。3.2 人为删失必须配权重否则等于自己给自己挖坑克隆造出来的人为删失不是随机的。谁会在某个策略下活下来即持续符合该策略跟他的基线特征和时变状态强相关。病情恶化的人更可能换药于是他们的原策略克隆被删失而留下来的人系统性地更稳定。直接对剩下的样本做分析等于在比较两组精英偏倚比不克隆还大。解决办法是给每个克隆按时间累积一个逆概率权重。直观理解如果一个克隆存活到 t 时刻的概率只有 0.2那它一个人就要代表五个像它这样的人权重就是 5。权重随时间累积公式大致是$$W_i(t) \prod_{k0}^{t} \frac{1}{P(C_i(k)0 \mid \bar{L}_i(k-1), A_i)}$$其中 $C_i(k)0$ 表示在 k 时刻仍未被删失$\bar{L}_i(k-1)$ 是到 k-1 时刻的历史协变量。分母用模型估计通常是 pooled logistic每个时间片拟合一次分子的选择决定了权重是原始的 1/P 还是稳定化的。稳定化权重是必做的。分子换成边际概率 $P(C_i(k)0 \mid C_i(k-1)0, A_i)$也就是不挂协变量的版本。这样权重的均值会被拉回 1 附近方差显著降低而点估计不变。我几乎没见过不稳定的权重能跑出可信结果的场景除非样本量特别大且正値性极好。实际写代码时权重通常分两块算再相乘策略依从权重处理换药、停药导致的人为删失和失访权重处理真实失访和行政截尾。两块分开算、分开诊断出了问题好定位。下面是一段简化示意重点是结构而不是具体变量library(speedglm) # 把数据展开成人-时间格式后逐时间片拟合删失模型 fit_censor - speedglm( censor_event ~ strategy age_at_t0 sex comorbidity_score time_since_t0 factor(time_slice), family binomial(), data dt_person_time ) dt_person_time[, p_uncensored : 1 - predict(fit_censor, type response)] # 稳定化分子用只含 strategy 与时间的简化模型 fit_num - speedglm( censor_event ~ strategy factor(time_slice), family binomial(), data dt_person_time ) dt_person_time[, p_num : 1 - predict(fit_num, type response)] # 按个体累积相乘得到稳定化权重 dt_person_time[, sw : cumprod(p_num / p_uncensored), by clone_id]有个实操细节必须提cumprod那一步一定要按clone_id分组不是按id。同一个个体的两个克隆走的是两条不同的路径权重各自累乘混在一起算会直接把结构搞乱。我第一次写的时候就在这里栽过结果权重分布完全对不上排查了大半天才发现是by写错了。3.3 权重诊断截尾位置、有效样本量和你该有的心理预期权重算出来之后第一步不是跑结局模型是看分布。我固定会看的几个数字稳定化权重的均值应该接近 1偏离太多说明分子模型有问题、最大值、99 分位数、以及有效样本量。有效样本量ESS的公式是 $(\sum w)^2 / \sum w^2$。它的含义很直观虽然你有 N 个克隆但由于权重不均实际起作用的样本量大约相当于 ESS 这么多。这个数字比 N 小得多的时候你就要警惕了。ESS 相对 N 的比例我的判断 50%权重比较健康可以正常分析20% - 50%需要检查极端权重考虑截尾5% - 20%结论可信度下降亚组分析基本别碰 5%停下来回头检查模型设定或策略定义截尾是最常用的补救手段。把权重在 1% 和 99% 分位或者固定阈值比如 10处截断用截断后的值替代原值。代价是引入少量偏倚换来方差大幅下降。截尾位置要在方案里预先写定不能跑完看结果不显著就往下调那样等于变相挑结果。我的习惯是同时报截尾和未截尾两版如果差别在可接受范围内说明结果对截尾不敏感。还有一点值得说权重偏大往往不是权重方法的问题是策略定义或协变量集的问题。比如你把一个用药后才会变化的变量放进了删失模型的分母那模型本身就在用未来的信息解释现在的行为权重结构必然扭曲。时变协变量要不要进模型判断标准是它是否同时影响后续是否继续该策略和结局。这个判断需要领域知识不能全交给统计软件。4. 估计量与置信区间从权重到最终那一行数字4.1 结局类型决定估计路径选错了会白干加权完之后用什么估计量取决于结局类型。这三条路我在不同项目里都走过各自的坑不太一样。时间到事件结局最常用。两条路一是加权后的乘积极限估计直接对每个策略组算加权生存曲线二是 pooled logistic 回归把数据展开成人-时间格式每个时间片一行拟合在 t 时刻之前未发生事件的前提下t 时刻发生事件的概率然后从模型里推出累积 Incidence。后者更灵活可以放时变协变量、可以直接算风险比和风险差也能顺带把置信区间算出来。代价是时间片划分要合理切太细样本稀、切太粗丢信息我一般按月或者按季度切看随访总长度。# person-time 格式下拟合 pooled logistic fit_out - speedglm( event ~ strategy factor(time_slice) age_at_t0 sex comorbidity_score, family quasibinomial(), data dt_person_time_w, weights sw # 稳定化权重进来用稳健方差 ) # 从模型系数推累积风险通常需要标准化predict 后取平均quasibinomial而不是binomial是为了避免权重不是整数时 R 弹出的警告。虽然点估计一样但后续算方差时的处理不同这个细节很多人会忽略。二分类结局比如随访期内是否发生某个不良事件直接算加权风险差和风险比就行每个策略组里加权事件数除以加权总人数。简单但要注意分母是权重和不是行数。用行数去除是新手最常犯的错误之一结果会明显偏离。连续结局比如某个指标的变化量用加权均值差。这里要注意的是如果结局在随访期内有多次测量取最后一次这种做法会引入删失偏倚活得更久的人有更多测量机会稳妥的做法是预先定义一个固定的测量时点窗口。4.2 方差估计软件默认输出基本都不能直接信权重是个双刃剑它修正了偏倚但让独立同分布这个前提彻底失效。软件看到权重通常当成频率权重处理算出来的标准误是偏小的。这在 TTE 里是系统性错误不是偶发。两个正经办法。一是稳健方差sandwich estimator很多包在拟合时加个选项就能出。二是bootstrap但有个致命细节必须按个体重抽样不能按行或按克隆重抽样。因为同一个个体的多个克隆之间高度相关按行抽样会把相关性打散得到的区间偏窄。正确的做法是抽id把抽中的这个人的所有克隆、所有时间片整块带进来。set.seed(20240101) ids - unique(dt$id) boot_est - replicate(500, { samp - sample(ids, length(ids), replace TRUE) d - dt[id %in% samp] # 整块抽取克隆和时间片一起带走 refit_pipeline(d) # 重跑权重 结局模型 }) # 取 2.5% 和 97.5% 分位作为 95% 区间 quantile(boot_est, c(0.025, 0.975))重跑 500 次完整的权重加结局模型在几十万行的数据上可能要跑几个小时。实际操作里我会先用 50 次看看分布形状对不对确认没有明显异常再放到 500 次。如果时间实在紧张至少也要 200 次。少于 100 次的 bootstrap 区间不稳定报上去自己都不放心。4.3 把可能错在哪量化出来阴性对照与 E-value任何观察性因果分析都躲不开未测量混杂这个质疑。与其在讨论部分写一段空话不如用两个工具把它量化。阴性对照分两种。阴性对照结局是找一个已知不受干预影响的结局比如与药物机制无关的某类意外伤害如果分析显示有效应那基本可以判定存在残余混杂或者别的系统性偏差。阴性对照暴露是找一个已知不影响目标结局的暴露用同样的流程跑一遍如果跑出效应同样说明流程有问题。这两个检查做起来不复杂但说服力很强我在正式分析里都会加。E-value给的是一个更直接的答案需要多强的未测量混杂才能把观察到的效应完全解释掉。公式是 $E RR \sqrt{RR \times (RR - 1)}$风险比小于 1 时先取倒数。算出来的数字越大结论越稳健。比如 RR 1.5 对应 E-value 约 2.37意思是需要一个与处理和结局的关联强度都达到 2.37 倍的未测量混杂才能把效应归零。把这个数字放在讨论里比写可能存在残余混杂有用得多。e_value - function(rr) { if (rr 1) rr - 1 / rr rr sqrt(rr * (rr - 1)) } e_value(1.5) # 2.37 e_value(2.0) # 3.41还有一个更进阶的做法是定量偏倚分析假设未测量混杂与已测量混杂的关联结构类似推算出校正后的效应范围。这套方法在方案里可以作为预设的敏感性分析写进去做不做取决于项目的重要程度但只要写了审稿人那关会好过很多。5. 落地时最容易翻车的几个细节5.1 想清楚要的是分配效果还是执行效果这是开工前必须回答的问题而且答案会改变整个分析流程。分配效果对应 ITT回答的是如果给一群人推荐 A 方案、给另一群推荐 B 方案实际执行中允许偏离最终结局差多少。执行效果对应 per-protocol回答的是如果所有人都严格按照方案执行结局差多少。克隆-删失-加权这套机制天然适合估执行效果因为你是在人为地把偏离策略的路径删失掉。如果你想要的是分配效果那就不该在换药时删失而是应该保留原策略标签继续随访——但这样又回到了分组依据是起始方案的老路需要额外的处理。很多文章在这两个目标之间含糊其辞读起来很别扭写的时候一定要在方案里点明。5.2 删失事件和结局事件的边界要提前划清哪些情况算删失、哪些算结局、哪些算竞争事件这三者的划分会显著影响结果而且没有放之四海皆准的答案。我的判断原则是如果这个事件发生后再观察原结局就没有意义了那它应该被当作竞争风险处理如果它只是让观察中断但理论上还能继续观察那是删失。举个具体的随访期内死亡。如果结局是某指标改善人死了当然不可能改善这时候死亡是竞争事件简单删失会高估改善率。如果结局是全因死亡那死亡就是结局本身。这两种处理方式在代码上就一行之差但解释完全不同方案里必须写清楚。还有一个隐蔽的坑行政截尾。数据截止到某一天之后的人为删失不是随机的越晚入组的人随访越短。加权模型里要把日历时间放进去否则随访长度和入组时间相关会引入偏倚。5.3 权重轰炸之后别被显著两个字骗了跑完看到 p 0.05先别高兴。TTE 里有一类假阳性特别典型权重把有效样本量压到只剩几百但软件输出的标准误还是按几十万行算的于是区间窄得离谱什么都显著。这种情况在亚组分析里尤其常见——整体 ESS 还有几千切到某个亚组只剩一两百结论完全不可靠。我的做法是任何亚组分析都单独报 ESSESS 低于某个阈值我一般用 200的亚组直接不报或者只做描述性呈现。整体结论对亚组是否一致用交互项检验而不是分开展示一堆单独显著的估计值。分开报容易给人处处有效的错觉实际上是多重比较加权重稀释的双重问题。5.4 报告清单审稿人和读者一定会盯的四个点最后整理一份我会固定检查的清单不管投哪里这四项都要能一句话答上来。检查项必须写清楚的内容时间零点具体到哪一天入组/分配/随访三者如何对齐策略定义起始条件、允许的偏离、换药如何处理权重构造模型里放了什么、稳定化分子是什么、截尾在哪里敏感性分析做了哪几项、结果是否一致、不一致时怎么解释前两项决定研究的可信度后两项决定审稿人会不会追问。我个人的经验是把这张表填满返修意见里关于方法不透明的部分基本就消失了。还有一个容易被忽略的补充材料代码和数据字典。TTE 的流程链条长任何一个环节换个写法结果都可能变。把代码整理干净放出去别人能复现你的结论才站得住。我现在的习惯是主分析脚本从原始数据读入开始一路到最后一张图和表中间不手工干预这样出了问题能从头查。如果一定要说一句我自己最深的体会TTE 的难点不在统计在把研究问题说清楚到能写成代码的程度。大多数做不下去的项目卡住的都是时间零点到底是哪天策略的边界在哪里这类看起来不像技术问题的技术问题。把这两件事想明白剩下的加权和建模反而是最省心的部分。