临床试验中重复测量数据的多重填补：SAS实现与应用

1. 临床试验中的重复测量数据缺失问题

在药物临床试验中，重复测量数据是最常见的数据类型之一。想象一下，我们需要定期测量患者的血压、血糖或其他关键指标，这些数据通常会按照预定的时间点（比如基线期、治疗第7天、第14天等）进行收集。但现实情况是，患者可能因为各种原因错过某些访视，导致数据出现缺失。

这种缺失可不是简单的数据丢失，它可能直接影响研究结论的可靠性。举个例子，如果某个治疗组的患者因为副作用提前退出研究，他们的缺失数据就可能带有特定模式。这时候，传统的统计方法（比如直接删除缺失数据）就会产生偏差。

我处理过的一个实际案例中，研究团队发现使用简单删除法得出的治疗效果比实际情况乐观了将近15%。这就是为什么我们需要更专业的缺失数据处理方法——多重填补（Multiple Imputation）。这种方法不是简单地猜测缺失值，而是通过建立概率模型，生成多个可能的填补结果，最后综合这些结果得出更可靠的结论。

2. 数据准备与模拟

2.1 创建完整的模拟数据集

在开始处理缺失数据前，我们需要一个基准数据集。使用SAS模拟数据是个好方法，因为我们可以完全控制数据的特征。下面这段代码模拟了一个典型的临床试验数据集：

PROC FORMAT; VALUE vis 1 = 'BL' 7 = 'DAY 7' 14= 'EOT' 28= 'FU D28' 42= 'FU D42' 98= 'FU D98'; RUN; DATA dd1; LENGTH subjid trtp $7. paramcd $4. param $17.; ARRAY visits (6) _TEMPORARY_ (1 7 14 28 42 98); CALL STREAMINIT(1977); DO id = 1 TO 200; subjid='101-'||PUT(id, Z3.); IF RAND('Bernoulli', 0.5) = 1 THEN sex='M'; ELSE sex='F'; IF RAND('Bernoulli', 0.5) = 1 THEN trtp='Active'; ELSE trtp='Placebo'; age=ROUND(RAND('Normal', 50, 15)); paramcd='MIDI'; param='Midichlorians (n)'; DO i=1 TO DIM(visits); avisitn=visits(i); avisit=PUT(avisitn, VIS.); IF trtp='Active' THEN DO; IF visits(i)=1 then aval=ROUND(RAND('Normal', 10000, 100)); ELSE aval=ROUND(RAND('Normal', 10000, 100)+visits(i)); END; ELSE DO; aval=ROUND(RAND('Normal', 10000, 100)); END; OUTPUT; END; END; DROP i id; RUN;

这段代码生成了200名患者的数据，包括：

• 患者ID（subjid）
• 治疗组别（trtp，分为Active和Placebo）
• 性别（sex）
• 年龄（age）
• 6个时间点的测量值（aval）

2.2 模拟缺失数据

现实中，数据很少是完整的。我们可以模拟不同类型的缺失模式：

DATA adeff; SET dd1; BY subjid; CALL STREAMINIT(1980); RETAIN base .; IF first.subjid THEN base=.; IF avisit='BL' THEN base= aval; IF avisit ne 'BL' THEN DO; IF RAND('Bernoulli', 0.08) THEN DELETE; END; RUN;

这里我们设置了8%的随机缺失率（非基线访视）。在实际分析中，了解缺失模式至关重要。常见的缺失机制有三种：

1. 完全随机缺失（MCAR）：缺失与任何变量无关
2. 随机缺失（MAR）：缺失与观察到的数据有关
3. 非随机缺失（MNAR）：缺失与未观察到的数据有关

3. 使用MMRM模型分析

3.1 混合效应模型重复测量（MMRM）

当数据是随机缺失（MAR）时，MMRM是个不错的选择。它直接利用现有数据建模，不需要填补缺失值。这种方法特别适合连续型结局变量。

data ana_p; set adeff; where avisitn in (7,14,28,42,98); chg=aval-base; run; proc mixed data=ana_p(where=( chg ne .)); class trtp(ref="Placebo") avisitn sex subjid; model chg=base sex age trtp avisitn trtp*avisitn /htype = 3 ddfm=kenwardroger; repeated avisitn / subject = subjid type=un; lsmeans avisitn*trtp / cl diff; ods output lsmeans=pls; run;

这段代码做了以下几件事：

1. 计算每个时间点相对于基线的变化值（chg）
2. 建立混合效应模型，考虑治疗组、性别、年龄等因素
3. 使用非结构化协方差矩阵（type=un）处理重复测量间的相关性
4. 输出最小二乘均值估计结果

3.2 MMRM的优缺点

优点：

• 直接利用现有数据，不需要填补
• 对随机缺失数据提供无偏估计
• 可以灵活处理不同协方差结构

缺点：

• 对非随机缺失（MNAR）可能产生偏差
• 结果解释相对复杂
• 不适合极端高缺失率的情况

4. 多重填补的实现步骤

4.1 数据格式转换

多重填补前，通常需要将数据从长格式转为宽格式：

PROC TRANSPOSE DATA=adeff OUT=onepersub PREFIX=MIDI; BY subjid age sex trtp base paramcd param; ID avisitn; VAR aval; RUN; proc means data=onepersub nmiss noprint; var MIDI1 MIDI7 MIDI14 MIDI28 MIDI42 MIDI98; output out=miss(drop=_type_ _freq_) nmiss=; run;

转换后的数据结构更便于填补，每个患者一行，各时间点的测量值作为不同变量。

4.2 多重填补过程

多重填补通常包括三个步骤：

1. 填补：生成多个完整数据集
2. 分析：对每个数据集分别分析
3. 合并：综合所有分析结果

4.2.1 判断缺失模式

PROC MI DATA=onepersub NIMPUTE=0; CLASS sex trtp; FCS; VAR sex trtp age base midi7 midi14 midi28 midi42 midi98; RUN;

这个过程帮助我们了解缺失是单调的还是任意的，以及各变量间的缺失关系。

4.2.2 执行多重填补

对于非随机缺失（MNAR），我们可以使用delta-adjusted方法：

proc mi data=onepersub out=imputed seed=125 nimpute=50; class sex trtp; var sex trtp age base midi7 midi14 midi28 midi42 midi98; fcs reg; mnar adjust(midi7 / shift=1 adjustobs=(trtp="Active")) adjust(midi14 / shift=1 adjustobs=(trtp="Active")) adjust(midi28 / shift=1 adjustobs=(trtp="Active")) adjust(midi42 / shift=1 adjustobs=(trtp="Active")) adjust(midi98 / shift=1 adjustobs=(trtp="Active")); run;

这里我们：

• 生成了50个填补数据集（nimpute=50）
• 对Active组的后期访视测量值进行了delta调整（shift=1）
• 使用FCS回归方法进行填补

4.2.3 分析填补后的数据

对每个填补数据集进行分析：

data imputed; set imputed; chg=midi98-base; run; proc mixed data=imputed; by _imputation_; class sex trtp(ref="Placebo"); model chg =base sex age trtp/htype=3 ddfm=kr; lsmeans trtp /diff alpha=0.05; ods output tests3 = etest3 diffs=ediffs lsmeans=els; run;

4.2.4 合并分析结果

proc sort data=els;by trtp;run; proc mianalyze data=els; by trtp; modeleffects estimate; stderr stderr; run; proc mianalyze data=ediffs; modeleffects estimate; stderr stderr; run;

这个过程使用Rubin规则合并估计值和标准误，给出最终的综合结果。

5. 实际应用中的注意事项

5.1 填补方法的选择

根据我的经验，选择填补方法时要考虑：

1. 缺失机制（MCAR/MAR/MNAR）
2. 变量类型（连续/分类）
3. 缺失模式（单调/任意）
4. 样本量和缺失比例

常用的SAS填补方法包括：

• 回归法（REG）
• 预测均值匹配（PMM）
• 马尔可夫链蒙特卡洛（MCMC）
• 完全条件设定（FCS）

5.2 填补质量的评估

填补后，应该检查：

1. 填补值的分布是否合理
2. 原始数据与填补数据的关系
3. 不同填补数据集间的变异性

可以使用PROC MI的diagnostic选项进行评估：

proc mi data=onepersub nimpute=5 out=mi_out; var midi7 midi14 midi28 midi42 midi98; fcs plots=trace; run;

5.3 敏感性分析

多重填补结果应该与主分析方法（如MMRM）进行比较。如果结果差异较大，可能需要：

1. 检查缺失机制假设
2. 调整填补模型
3. 考虑不同的delta调整值

我曾遇到一个案例，当delta值从1增加到1.5时，治疗效果估计值变化了12%，这说明结果对MNAR假设相当敏感。

6. 代码优化与实用技巧

6.1 提高计算效率

多重填补可能很耗时，特别是当：

• 数据集大
• 填补次数多
• 模型复杂

几个优化建议：

1. 使用BY语句替代CLASS语句
2. 简化协方差结构
3. 合理设置nimpute值（通常20-100次足够）

6.2 处理大规模数据

对于大型临床试验数据：

1. 分阶段填补：先填补关键变量，再填补次要变量
2. 使用SAS/STAT的高性能程序
3. 考虑并行处理

6.3 常见错误排查

我总结的几个常见问题：

1. 分类变量未正确指定
2. 种子值设置不当导致结果不稳定
3. 缺失模式诊断不充分
4. 合并步骤忽略了BY组处理

一个实用的调试技巧是先用小样本（nimpute=3）测试代码，确认无误后再运行完整分析。