从矩阵基础到监管提交

1.1 什么是纵向数据？

在临床试验中，我们通常不会只在终点测量一次结局指标，而是在多个访视（visit）对同一受试者进行重复测量。例如，一项疼痛研究可能在基线（Baseline）后第 1、2、4、6、8、12 周分别评估疼痛评分（NRS）。这种同一受试者在多个时间点的重复测量就是纵向数据（longitudinal data）。

1.2 纵向数据的三大难题

方法对比

4.1 模型形式

设受试者 $i$ 在访视 $j$ ($j = 1, \ldots, k$) 的相对基线变化为 $Y_{ij}$：

$$Y_{ij} = \mu + \tau_{t(i)} + \gamma_j + (\tau\gamma)_{t(i)j} + \beta \cdot \text{base}_i + \varepsilon_{ij}, \quad \boldsymbol{\varepsilon}_i \sim N(\mathbf{0}, \boldsymbol{\Sigma})$$

4.2 混合模型 vs MMRM

「混合模型」(Mixed Model) 是指同时包含固定效应和随机效应的模型。MMRM 是一种特殊的混合模型 —— 它通常只有 R 侧（残差协方差），没有 G 侧（随机效应协方差）。

4.3 交互模型构建器

选择模型项，实时查看模型方程和 SAS 代码如何变化：

5.1 MODEL 和 REPEATED：一张表的两个方向

沿用一个虚构研究框架：受试者在 W1、W2、W4、W6、W8 和 W12 评估疼痛 NRS，主要分析因变量是 $\text{CHG} = \text{AVAL} - \text{BASE}$。

MODEL 语句描述的是：给定基线、治疗组、访视和分层因素，平均的 CHG 处在什么位置？

REPEATED 语句描述的是：同一个人的多个 CHG 如何一起变化？

5.2 常见结构对比

5.3 SUBJECT=USUBJID：先划定「同一个人」

在 REPEATED AVISITN / SUBJECT=USUBJID TYPE=UN; 中，SUBJECT=USUBJID 告诉 SAS：具有相同 USUBJID 的记录属于同一个重复测量簇。

程序员需要检查：

• USUBJID 是否真正唯一标识受试者
• 同一受试者是否被错误拆成多个 ID
• 是否把治疗组、中心或访视误写成 SUBJECT=
• 每位受试者每个分析访视是否最多一条目标记录
• REPEATED 中的访视变量是否与数据排序和分析访视一致

5.4 协方差结构交互实验室

选择不同结构、调整访视数和相关参数，实时观察协方差矩阵的变化：

6.1 三种缺失机制

6.2 MMRM 如何处理缺失：基于似然，无需插补

这句话包含四个要点，每一条都对应一个具体的建模或编程动作：

之所以成立，关键在于 $V$ 是按每位受试者实际观测到的访视取子块：第 $i$ 名受试者贡献的似然只涉及 $V_i$ 中对应的行与列，缺失访视对应的行列被直接划掉，而不必先猜一个值填进去。因此 MMRM 与「先插补再分析」是两条不同的技术路线——后者人为增加了信息量，会低估不确定性，通常导致标准误偏小、P 值偏小。

6.3 缺失机制交互演示

观察不同缺失机制下数据点的缺失模式差异：

6.4 MMRM 不会生成单次填补值

这句话容易在程序复核时被忽略。把"无填补"拆成三个可核查的命题：

三种"缺失"在数据结构上不同，但都不会产生虚构的结局值：

从似然的角度看（推导见 S23）：每名受试者只按自己实际的观测访视提供信息，$V_i = Z_i \Sigma Z_i^\top + \sigma^2 I_{n_i}$ 中的 $n_i$ 因人而异。模型评估的是"在当前均值模型与协方差结构下，已观测到的结果与模型有多一致"，而不是"缺失值该填什么"。

程序员可以核验的六项：

6.5 为什么 MAR 主分析之外仍需敏感性分析

MMRM 在 MAR 下能给出无偏估计与有效推断，但这不等于 MNAR 风险已经消失。关键在于：

因此 ICH E9(R1) 要求：缺失数据的处理必须相对于具体的 estimand 来定义，并用预先设定的敏感性分析检验主分析对关键假设偏离是否稳健。主分析与敏感性分析面向同一 estimand，通过改变"缺失结局的假设"来观察结论是否稳定。

常见的参考组敏感性分析策略：

J2R 与 CR 的概念差异（假设某受试者在 W6 后退出，安慰剂组为参考组）：

6.6 MNAR 敏感性分析的生产实现：两段式 MI + 参考组填补

下面这段代码取自另一份真实 SAP（文中统称 STUDY-XX-202），是"MMRM 主分析 + MNAR 敏感性分析"配套写法中最典型的一种：先用 MAR 假设把间歇缺失补成单调缺失模式，再用参考组（control-based）FCS 在 MNAR 假设下补齐剩余的单调缺失。

先明确它挂在哪个 estimand 上 —— 这是程序员最容易搞混的地方：

SAP 中的示例程序（变量名与种子保留原样，研究号与药物代号已替换）：

proc mi data=datain_wide
        out=mi_mono
        nimpute=100
        seed=201101;
  var TRT STRAT BASE V2 V4 V8 V12;
  mcmc impute=monotone;
run;

proc sort data=mi_mono;
  by _imputation_;
run;

proc mi data=mi_mono
        out=mi_fcs
        nimpute=1
        seed=201102;
  by _imputation_;
  class TRT STRAT;
  var BASE STRAT TRT V2 V4 V8 V12;
  fcs
    reg( V2  = BASE STRAT TRT )
    reg( V4  = BASE STRAT TRT V2 )
    reg( V8  = BASE STRAT TRT V2 V4 )
    reg( V12 = BASE STRAT TRT V2 V4 V8 );
  mnar model ( V2 V4 V8 V12 / modelobs=(TRT='PBO') );
run;

随后对 mi_fcs 中每个 _imputation_ 运行与主分析完全相同的 MMRM，用 PROC MIANALYZE（Rubin's rules）合并 100 组 LSMean 差值与 SE。

逐行解读：

报告与呈现的要求：

• 主分析与敏感性分析并列呈现在同一张表中（点估计、SE、CI、P 值四列对齐），让读者直接比较；
• 结论稳健性按 SAP 预先定义的标准判断，例如"主分析与敏感性分析方向一致且结论不变"；
• 若两者结论不一致，不能选择性报告，必须在 CSR 中说明差异来源（缺失率、退出原因分布、插补模型）；
• 二分类终点常用 NRI（non-responder imputation）处理缺失，那是另一套机制，不要与连续终点的 MNAR 插补混为一谈。

7.1 REML vs ML

最大似然（ML）和限制性最大似然（REML）是两种估计协方差参数的方法：

SAS 中指定 REML：

PROC MIXED DATA=adqs METHOD=REML;
  CLASS TRTP AVISIT USUBJID;
  MODEL CHG = TRTP AVISIT TRTP*AVISIT BASE / DDFM=KR;
  REPEATED AVISIT / SUBJECT=USUBJID TYPE=UN R RCORR;
  LSMEANS TRTP*AVISIT / DIFF CL;
RUN;

7.2 自由度校正方法

左边是一份 ADaM BDS 长格式数据（3 名受试者 × 4 个访视，含 2 个缺失），右边是 MMRM 的完整代码。点「下一步」逐步推进：亮起的列就是这一步真正被读取的数据，曲线连到代码里对应的关键字，下方卡片解释这一步在数学上做了什么运算。

9.1 RANDOM vs REPEATED：G 侧与 R 侧

$$\mathbf{V}_i = \mathbf{Z}_i\mathbf{G}\mathbf{Z}_i' + \mathbf{R}_i \quad \text{（RANDOM 建 G，REPEATED 建 R）}$$

10.1 参数数量计算器

输入访视数，查看各结构需要估计的协方差参数数量：

10.2 何时选择何种结构？

选择下面的参数组合，查看对应的模型规格和预期输出：

11.1 每个输出参数是怎么算出来的

沙盘里切换的是「模型设置」，而输出表格里的每一个数字都有确定的计算路径。下表把 SAS 输出中的关键量、所用公式、输入数据与估计方法一一对应，便于逐项复核。

11.2 关键中间结果：用 W16 差值做一次完整复核

下面用 S22 那套真实输出（STUDY-XX-201 某亚组，CS 结构，DDFM=KR）把上表公式逐项算一遍。统计量数值未作改动，仅研究编号与治疗组代号已替换。你自己复核时，只要拿到 CovParms、LSMeans 与 Diffs 三张表就能重现以下全部数字。

Step 1 · 协方差参数 → $V$ 的结构

CS 结构下，$\hat V_i$ 的对角元是两个参数之和，非对角元为 CS 分量：

• $\operatorname{Var}$（对角）$= \widehat{CS}+\widehat{Residual}=73.2028+61.7798=134.9826$
• $\operatorname{Cov}$（非对角）$= \widehat{CS}=73.2028$
• 隐含的组内相关 $= 73.2028\,/\,134.9826 = 0.5423$

Step 2 · LSMean 与治疗差值

• 研究药物 W16 LSMean $= -43.2597$（SE 8.9220，DF 48.7）
• 安慰剂 W16 LSMean $= -18.4896$（SE 9.0975，DF 53.1）
• 差值（研究药物 − 安慰剂）$= (-43.2597)-(-18.4896) = \mathbf{-24.7701}$，Diffs 表报 $\mathbf{-24.7700}$（末位舍入），两者一致

注：早期版本此处引用过另一组 LSMean/SE 数值，来源为对纯文本输出的错误解析（列错位）。本版已按原始 sashtml.htm 的 Least Squares Means 与 Differences of Least Squares Means 两表逐格核对更正。

Step 3 · SE(差值)：协方差项不能丢

• $\operatorname{Var}_1 = 8.9220^2 = 79.6021$；$\operatorname{Var}_2 = 9.0975^2 = 82.7645$
• $\operatorname{Var}(\text{diff}) = 12.7445^2 = 162.4223$
• $2\operatorname{Cov}_{12} = 79.6021+82.7645-162.4223 = -0.0557 \Rightarrow \operatorname{Cov}_{12}=-0.0278$
• $\operatorname{Corr}_{12} = -0.0278\,/\,(8.9220\times9.0975) = -0.0003$

Step 4 · t 值、置信区间与 P 值

• $t = -24.7700\,/\,12.7445 = -1.9436$（SAS 显示 −1.94）
• $t_{0.975,\,51.7} \approx 2.0071 \Rightarrow$ 边际误差 $= 2.0071\times12.7445 = 25.580$
• 95% CI $= (-24.7700-25.580,\ -24.7700+25.580) = (-50.35,\ 0.81)$；SAS 输出 $(-50.3478,\ 0.8077)$，差异仅来自分位数的四舍五入
• $p = 2\big[1-T_{51.7}(1.9436)\big] = 0.0574$

Step 5 · 拟合统计量的 $n$ 口径

已知 $N=80$，$\operatorname{rank}(X)=11$，$d=2$（CS），受试者数 $=36$，−2RLL $=545.711$：

• AIC $= 545.711 + 2\times2 = 549.711$（SAS 549.7）✓
• AICC：$n = N-\operatorname{rank}(X)=80-11=69$，$\Rightarrow 545.711 + \dfrac{2\times2\times69}{69-2-1} = 545.711+4.1818 = 549.893$（SAS 549.9）✓
• BIC：$545.711 + 2\times\ln 36 = 545.711 + 7.1670 = 552.878$（SAS 552.9）✓

三者的 $n$ 不同：AIC 不用 $n$；AICC 用 $N-\operatorname{rank}(X)$；BIC 用受试者数。拿错口径就复现不出 SAS 的数字——这是自查时最容易卡住的地方。

11.3 参数复核清单

1. 先对结构再看数字。 CovParms 的参数个数是否等于该 TYPE= 的理论个数？出现 0、负值或贴边界，先回到 S24 排查。
2. 用 E 选项打印 $\ell$ 向量。 确认 LSMean / 差值确实用了你以为的那一组系数，尤其是有协变量或交互项时。
3. 差值必须来自 Diffs 表。 不要用取整后的两个 LSMean 相减；同样的道理，CI 与 P 值也必须用模型直接给出的 SE 与 DF 重算。
4. 检查 SE 与 Var 的自洽性。 用 Step 3 的式子反算 $\operatorname{Cov}_{12}$，若为负得不合常理，多半是 $\ell$ 向量或模型设定有问题。
5. DF 非整数不是错。 KR 的 $m$ 本来就是实数；但 DF 明显小于受试者数时要警觉单元格稀疏。
6. 用 AIC / AICC / BIC 反查 $n$ 口径。 三者能对上 SAS，说明 −2RLL 与 $d$ 都取对了。
7. 交叉验证。 条件允许时用另一套实现（如 R mmrm 或本指南的自研 REML 引擎）重跑同一模型，比对到小数点后 2–3 位。

12.1 SAS PROC MIXED 主分析模板

/*================================================================
  MMRM 主分析程序模板
  数据集：ADaM BDS 长格式 (ADQS)
  因变量：CHG = AVAL - BASE
  固定效应：TRTP, AVISIT, TRTP*AVISIT, BASE
  协方差：UN (非结构化)
  自由度：Kenward-Roger
  ================================================================*/

* 1. 数据准备：排序确保受试者内记录连续且按访视顺序;
PROC SORT DATA=ADQS OUT=ADQS_SORTED;
  BY USUBJID AVISITN;
RUN;

* 2. 受试者-访视唯一性检查;
PROC SORT DATA=ADQS_SORTED NODUPKEY;
  BY USUBJID AVISITN;
RUN;

* 3. MMRM 主分析;
PROC MIXED DATA=ADQS_SORTED METHOD=REML;
  CLASS TRTP AVISIT USUBJID;
  MODEL CHG = TRTP AVISIT TRTP*AVISIT BASE
        / DDFM=KR SOLUTION CL;
  REPEATED AVISIT /
    SUBJECT=USUBJID
    TYPE=UN
    R RCORR;
  LSMEANS TRTP*AVISIT / DIFF CL ADJUST=NONE;
  ODS OUTPUT
    LSMeans=LSM_OUT
    Diffs=DIF_OUT
    CovParms=COV_OUT
    FitStatistics=FIT_OUT;
  TITLE "MMRM Primary Analysis - UN Covariance";
RUN;
QUIT;

* 4. 探索性协方差结构比较;
%MACRO FIT_COV(TYPE=, LABEL=);
PROC MIXED DATA=ADQS_SORTED METHOD=REML;
  CLASS TRTP AVISIT USUBJID;
  MODEL CHG = TRTP AVISIT TRTP*AVISIT BASE / DDFM=KR;
  REPEATED AVISIT / SUBJECT=USUBJID TYPE=&TYPE;
  ODS OUTPUT FitStatistics=FIT_&LABEL;
RUN;
QUIT;
%MEND;

%FIT_COV(TYPE=UN,    LABEL=UN);
%FIT_COV(TYPE=CS,    LABEL=CS);
%FIT_COV(TYPE=AR(1), LABEL=AR1);
%FIT_COV(TYPE=TOEP,  LABEL=TOEP);

12.2 R mmrm 等效代码

# MMRM 主分析 - R mmrm 包
library(mmrm)

# 拟合 MMRM 模型
fit <- mmrm(
  formula = CHG ~ TRTP * AVISIT + BASE,
  cov_structure = "us",          # 非结构化协方差
  data = adqs,
  method = "REML"
)

# 查看模型摘要
summary(fit)

# 分访视治疗组对比（LS Mean Differences）
lsmeans <- emmeans(fit, specs = pairwise ~ TRTP | AVISIT)
print(lsmeans)

# 探索性协方差结构比较
fit_cs   <- update(fit, cov_structure = "cs")
fit_ar1  <- update(fit, cov_structure = "ar1")
fit_toep <- update(fit, cov_structure = "toep")

# 比较 AIC/BIC
AIC(fit, fit_cs, fit_ar1, fit_toep)
BIC(fit, fit_cs, fit_ar1, fit_toep)

13.1 关键输出表

13.2 输出 → TLF 映射

14.1 十二个常见误区

14.2 自检清单

勾选以下项目确认你的 MMRM 分析设置正确：

15.1 估计目标的五个属性

15.2 五种伴发事件处理策略

16.1 监管机构对比

16.2 敏感性分析

ICH E9(R1) 强调：主要分析（MMRM under MAR）应辅以敏感性分析，评估结论对 MAR 假设的稳健性。常见的敏感性分析方法包括：

• Pattern Mixture Models：对不同缺失模式分别建模
• Selection Models：联合建模结局和缺失机制
• Tipping Point Analysis：探索 MNAR 程度达到多少时结论改变
• Multiple Imputation with Delta Adjustment：在 MI 框架下引入 MNAR 偏移

监管文件

• ICH E9(R1). Addendum on Estimands and Sensitivity Analysis in Clinical Trials. 2019.
• FDA. Guidance for Industry: The Use of Multiple Imputation in Clinical Trials with Missing Data. 2010.
• EMA. Reflection Paper on Missing Data in Confirmatory Clinical Trials. 2010.
• FDA. Guidance for Industry: E9 Statistical Principles for Clinical Trials. 1998.

统计方法

• Mallinckrodt CH, et al. Recommendations for the primary analysis of continuous endpoints in longitudinal clinical trials. Drug Inf J. 2008;42(4):303-319.
• Mallinckrodt CH, et al. Assessing and interpreting treatment effects in longitudinal clinical trials with missing data. Biol Psychiatry. 2008;63(8):772-778.
• Kenward MG, Roger JH. Small sample inference for fixed effects from restricted maximum likelihood. Biometrics. 1997;53(3):983-997.
• Littell RC, et al. SAS for Mixed Models. 2nd ed. SAS Institute; 2006.
• Diggle PJ, et al. Analysis of Longitudinal Data. 2nd ed. Oxford University Press; 2002.
• Rubin DB. Multiple Imputation for Nonresponse in Surveys. Wiley; 1987.

软件文档

• SAS/STAT User's Guide: The MIXED Procedure. SAS Institute.
• Kenward MG, Roger JH. The mmrm package for R. GitHub: openpharma/mmrm.
• R Core Team. R: A Language and Environment for Statistical Computing. R Foundation; 2024.

社区资源

• 微信公众号「生物统计与统计编程」MMRM 系列文章
• PHUSE / CDISC MMRM 最佳实践指南
• Stack Overflow: [sas] [mixed-models] [mmrm] 标签

SAS 代码速查

PROC MIXED DATA=adqs METHOD=REML;
  CLASS TRTP AVISIT USUBJID;
  MODEL CHG = TRTP AVISIT TRTP*AVISIT BASE / DDFM=KR;
  REPEATED AVISIT / SUBJECT=USUBJID TYPE=UN R RCORR;
  LSMEANS TRTP*AVISIT / DIFF CL;
RUN;

R 代码速查

library(mmrm)
fit <- mmrm(CHG ~ TRTP * AVISIT + BASE,
            cov_structure = "us", data = adqs)
summary(fit)
emmeans(fit, pairwise ~ TRTP | AVISIT)

20.1 先分清三层数据的职责

ODS 结果不宜直接改造成显示字符串。更稳妥的做法是分成三层，每层只承担一种责任：

20.2 ODS 表名与触发语句对照

SAS 官方文档的 ODS Table Names 表规定了"表名 ← 触发它的语句/选项"这一一对应关系。写 ODS OUTPUT 之前先确认触发条件，否则会得到空数据集或根本不生成：

20.3 LSMeans 与 Diffs 的字段

LSMeans 每行对应一个「治疗组 × 访视」单元的调整后均值：

Diffs 每行是一条两两比较，端点字段带下划线：

20.4 LSMESTIMATE 语句：直接对 LSMean 做线性组合

当 SAP 要求的估计量不是现成的某一条两两比较时（例如"两个剂量组合并后与安慰剂比较"、"某访视上若干单元格的平均"），用 LSMESTIMATE 直接写 LSMean 的系数最直观——它作用的单位是 LSMean，而不是固定效应参数 β。

proc mixed data=adeff_ana method=reml cl;
    class usubjid trt01pn(ref="1") avisitn;
    model chg = base trt01pn avisitn trt01pn*avisitn / ddfm=kr solution;
    repeated avisitn / subject=usubjid type=un;

    /* LSMean 的水平顺序：avisitn 在外层、trt01pn 在内层，
       即 (W2,T1)(W2,T4)(W4,T1)(W4,T4)... 共 10 个单元 */
    lsmestimate trt01pn*avisitn
        'W16: 研究药物 - 安慰剂'  0 0 0 0 0 0 0 0 -1 1,
        'W12 与 W16 差值之差'     0 0 0 0 0 0 -1 1 1 -1
        / cl divisor=1,1 joint e;
    ods output LSMEstimates=lsm_est Coef=lsm_coef;
run;

输出数据集 LSMEstimates 的常用字段：Effect（对应的 LSMEANS 效应）、Label（你写的行标签）、Estimate、StdErr、DF、tValue、Probt、Alpha、Lower、Upper。加了 JOINT 时还会额外输出联合检验（F 值与 NumDF/DenDF）。

20.5 SLICE：选项与语句的区别

这是一个高频混淆点，先给结论：

/* 写法一：LSMEANS 的 SLICE= 选项 */
lsmeans trt01pn*avisitn / slice=avisitn cl;
ods output Slices=slc_a LSMeans=lsm;

/* 写法二：SLICE 语句（可额外要差值） */
slice trt01pn*avisitn / sliceby=avisitn diff cl adjust=tukey;
ods output Slices=slc_b;

两者都会产出 ODS 表名 Slices，输出标题为 Tests of Effect Slices。常用字段：

20.6 STORE 语句：把模型存下来，交给 PROC PLM

STORE 不产生任何输出表。它把分析上下文与结果写入一个二进制的 item store，之后由 PROC PLM 读取并做后处理——不需要重新拟合模型。对于跑一次要几十分钟的 UN 结构，这是最实在的省时手段。

proc mixed data=adeff_ana method=reml cl;
    class usubjid trt01pn(ref="1") avisitn;
    model chg = base trt01pn avisitn trt01pn*avisitn / ddfm=kr solution;
    repeated avisitn / subject=usubjid type=un;
    lsmeans trt01pn*avisitn / cl;
    store mmrm_store / label='W16 主分析 MMRM（UN, DDFM=KR）';
run;

/* 后续任意次后处理，都不再重跑模型 */
proc plm restore=mmrm_store;
    show cov parms;                       /* 查看 item store 内容 */
    lsmeans trt01pn*avisitn / diff cl;    /* 继续用 KR 自由度 */
    lsmestimate trt01pn*avisitn 'W16 diff' 0 0 0 0 0 0 0 0 -1 1 / cl;
    slice trt01pn*avisitn / sliceby=avisitn diff;
run;

21.1 DIFFs 的规模：为什么是 45 行

LSMEANS effect / DIFF 请求的是该效应全部水平组合之间、无重复的两两比较。若交互效应有 k 个单元格，Diffs 的行数就是组合数：

$\text{行数} = \binom{k}{2} = \dfrac{k(k-1)}{2}$

以本指南贯穿使用的真实 MMRM 输出为例：2 个治疗组 × 5 个访视 = 10 个单元格，于是 Diffs 有

$\binom{10}{2} = \dfrac{10\times 9}{2} = 45\ \text{行}$

而这 45 行里，访视内的治疗组比较只有 5 行（每个访视 1 条）；其余 40 行是跨访视比较，例如「研究药物 W2 的 LSMean 减去安慰剂 W16 的 LSMean」，临床上通常没有解释意义。

21.2 差值方向由两端决定，反向时 CI 必须成套翻转

Estimate = 第一端 − 第二端。程序不应假定 ODS 一定按某个方向输出，稳妥做法是同时接受两个方向再统一：

data w16_diffs;
    set mmrm_diffs;
    where avisitn=16 and _avisitn=16          /* 两端都锁在 W16 */
      and (  (trt01pn in (4)   and _trt01pn=1)      /* 研究药物 - 安慰剂 */
          or (_trt01pn in (4)  and trt01pn=1) );    /* 安慰剂 - 研究药物（反向） */
    if trt01pn in (4) and _trt01pn=1 then do;
        trtn=trt01pn; diff=estimate; lcl=lower; ucl=upper;
    end;
    else do;
        trtn=_trt01pn; diff=-estimate; lcl=-upper; ucl=-lower;   /* 注意：上下限互换 */
    end;
    diff_se=stderr; diff_df=df; pvalue=probt;
run;

21.3 ESTIMATE / CONTRAST / LSMESTIMATE 三方对照

三者的本质区别在于作用对象：ESTIMATE 与 CONTRAST 作用于固定效应参数 $\boldsymbol\beta$，LSMESTIMATE 作用于最小二乘均值。

21.4 同一个 W16 差值的三种写法

目标估计量固定为：W16 时「研究药物 − 安慰剂」的 LSMean 差值。模型为 CHG = BASE TRT01PN AVISITN TRT01PN*AVISITN，trt01pn(ref="1")、AVISITN 参考水平为 16。

estimate 'W16 diff'
    trt01pn 1 -1
    trt01pn*avisitn
        0 0 0 0 0 0 0 0 1 -1
    / cl;

contrast 'W16 diff'
    trt01pn 1 -1
    trt01pn*avisitn
        0 0 0 0 0 0 0 0 1 -1;

lsmestimate trt01pn*avisitn
    'W16 diff' 0 0 0 0 0 0 0 0 -1 1
    / cl e;

21.5 各自适用场景

22.1 分析背景

22.2 实际运行的程序

/* 0. 清理上一轮结果：避免不收敛时残留旧数据集被误用 */
proc datasets lib=work nolist;
    delete lsm diffs;
quit;

/* 1. 先探 UN 是否收敛（只要 ConvergenceStatus，不要结果表） */
ods select none;
ods output ConvergenceStatus=conv_un;
proc mixed data=sub_ana(where=(chg>.)) method=reml;
    class usubjid trt01pn avisitn;
    model chg = base trt01pn avisitn trt01pn*avisitn / ddfm=kr residual;
    repeated avisitn / subject=usubjid type=un;
run;
ods select all;

/* 2. 若 UN 未收敛，遍历回退链，记录各结构 AIC */
%let methlist=TOEPH|ARH(1)|CS|CSH|SP(POW)(avisitn_n);
/*   …… 逐个 PROC MIXED，把 FitStatistics 中的 AIC 汇总到数据集 AIC …… */

/* 3. 按 AIC 最小选定结构，重跑并取结果 */
proc mixed data=sub_ana(where=(chg>.)) method=reml cl;
    class usubjid avisitn trt01pn(ref="1");
    model chg = base trt01pn avisitn trt01pn*avisitn / ddfm=kr solution;
    repeated avisitn / subject=usubjid type=&curr_type.;
    lsmeans trt01pn*avisitn / pdiff cl;
    ods output LSMeans=lsm Diffs=diffs ConvergenceStatus=conv;
run;

/* 4. 收敛判定：reason 必须正好是 "Convergence criteria met." */
proc sql noprint;
    select count(*) into :fail_cnt from conv
    where reason ne "Convergence criteria met.";
quit;

22.3 逐表解读

① Model Information —— QC 第一步

本例的第一条结论：主结构 UN 没有收敛，程序按 SAP 回退链落到了 CS。这个事实必须在 QC 记录和 CSR 脚注中体现——"用的是什么结构"本身就是需要申报的信息。

② Class Level Information —— 数清楚水平

③ Dimensions —— 19 与 11 的差别

④ Number of Observations —— 缺失了多少

36 名受试者 × 5 个访视 = 180 个理论观测点，实际只有 80 条进入分析，缺失 100 条（55.6%）。Read = Used 且 Not Used = 0 说明没有因协变量缺失被剔除的记录。这个缺失比例解释了为什么 UN（15 个协方差参数）在此亚组上无法收敛。

⑤ Iteration History —— 优化过程本身

• 目标函数单调下降：567.61 → 545.71，共下降 21.90。
• Criterion 是基于梯度的收敛判据（不是相邻两次目标函数的差），从 0.012 一路降到 $3\times10^{-8}$。
• Iteration 1 的 Evaluations = 2：这一步做了线搜索（步长折半），所以多算了一次似然。
• 判读要点：下降是否平稳。若出现反复震荡、或目标函数在最后几步仍在明显变化，即使打了"Convergence criteria met"也要警惕。

⑥ Covariance Parameter Estimates —— 方差成分

CS 结构的含义是 $\Sigma = \sigma_{cs}\,J + \sigma^2_e\,I$，因此：

• 各访视方差 $= 73.2028 + 61.7798 = 134.9826$，即 SD = 11.62
• 组内相关系数（ICC）$= \dfrac{73.2028}{134.9826} = \mathbf{0.5423}$
• 两条置信区间都不含 0，且 CS 的下限 24.14 远离 0，说明受试者内相关确实存在

⑦ Fit Statistics 与 Null Model LRT —— 三个指标用了三个不同的 n

Null Model Likelihood Ratio Test：DF = 1，$\chi^2$ = 21.90，$p < .0001$。它检验的是 $H_0:\ \sigma_{cs}=0$（即退化为独立残差）。1 个自由度、$\chi^2$ 显著，说明纳入受试者内相关是必要的。

⑧ Solution for Fixed Effects（节选）

⑨ Type 3 Tests of Fixed Effects

⑩ Least Squares Means（10 行 = 2 组 × 5 访视）

注意 Pr > |t| 这一列是该 LSMean 与 0 比较的 P 值（"该组的平均变化是否不为 0"），不是组间比较 P 值。想回答"研究药物是否优于安慰剂"，必须去 Diffs。

⑪ Differences of Least Squares Means —— 45 行里只取 5 行

如 21.1 所述，10 个单元格产生 $\binom{10}{2}=45$ 条比较。其中访视内的组间比较只有 5 条：

读结果的结论：W16 时研究药物相对安慰剂的 DAPSA 变化差值估计为 −24.77（负号方向有利于研究药物，因为 DAPSA 越低越好），95% CI 为 $(−50.35,\ 0.81)$ 跨过 0，双侧 $p=0.0574$。因此该亚组在 W16 未达到 0.05 水平的统计学显著性——但这是一个亚组分析，样本量小、缺失率高，本就不以支持确证性结论为目的。

22.4 独立复现：自研 REML 估计器的交叉验证

本机没有 SAS，为验证上述解读所依据的计算逻辑，我用 numpy/scipy 从零实现了一套 REML 估计器（解析 score 与 Fisher 信息、Fisher scoring 迭代、Kenward-Roger 校正），在相同设计骨架的模拟数据上重跑：36 名受试者、5 个访视、80 条观测、缺失率 55.6%，真实协方差参数就设为真实输出估计出的 CS=73.2028 / Residual=61.7798。

22.5 三个必须记住的陷阱

22.6 从统计量到临床解读：−24.77 到底意味着什么

把 LSMean 差值读成结论，需要三步，缺一步就会出错。下面用本例真实的访视别差值走一遍。

访视别差值：真实的完整轨迹

终点为 DAPSA 较基线变化，数值越低越好，因此负值表示研究药物组改善更多。数据来源：Differences of Least Squares Means，同访视 TRT01PN=4 减 TRT01PN=1。

轨迹能读出什么（以及不能读出什么）

W16 的 −24.77，临床结论该怎么写

23.1 模型与符号

MMRM 是没有随机效应的混合模型：$Z = 0$，$G$ 不出现，因此 $V = R$。对受试者 $i$：

$y_i = X_i\boldsymbol\beta + e_i,\qquad e_i \sim N(0,\ V_i)$

全部受试者堆叠后 $V = \mathrm{blockdiag}(V_1,\dots,V_n)$，$\mathrm{rank}(X)=p$，$N=\sum_i m_i$。

23.2 $V$ 的构造：缺失是怎么进入矩阵的

这是理解"MMRM 不填补缺失"的关键。设 $\Sigma$ 是 $T\times T$ 的完整访视协方差矩阵，受试者 $i$ 实际观测到的访视下标集合为 $\mathcal{O}_i$（如 $\{1,3,5\}$ 表示只有 W2、W8、W16 有数据），则

$V_i = \Sigma[\mathcal{O}_i,\ \mathcal{O}_i]$

即取出 $\Sigma$ 对应的行与列构成子矩阵。受试者 $i$ 贡献 $m_i=|\mathcal{O}_i|$ 行，不同受试者的 $m_i$ 可以不同。全程没有任何一步需要填补缺失值——这就是 MAR 下 MMRM 能直接使用不完整纵向数据的矩阵解释，也是它与 LOCF / MI 的根本区别。

23.3 从 ML 到 REML：把 $\boldsymbol\beta$ 积分掉

完整数据的对数似然为

$\ell_{ML}(\boldsymbol\beta,\theta) = -\dfrac{N}{2}\log 2\pi - \dfrac12\log|V| - \dfrac12 (y-X\boldsymbol\beta)'V^{-1}(y-X\boldsymbol\beta)$

ML 把 $\boldsymbol\beta$ 当作未知参数一起估计，导致方差参数被系统性低估（未扣除 $p$ 个自由度）。REML 的做法是把 $\boldsymbol\beta$ 从似然中"积分掉"，只基于 $N-p$ 个与 $\boldsymbol\beta$ 正交的误差对比做推断。

核心的代数分解（关键一步）：

$(y-X\boldsymbol\beta)'V^{-1}(y-X\boldsymbol\beta) = r'V^{-1}r + (\boldsymbol\beta-\hat{\boldsymbol\beta})'(X'V^{-1}X)(\boldsymbol\beta-\hat{\boldsymbol\beta})$

其中 $\hat{\boldsymbol\beta} = (X'V^{-1}X)^{-1}X'V^{-1}y$ 是 GLS 估计，$r = y - X\hat{\boldsymbol\beta}$。右侧第二项与 $\boldsymbol\beta$ 有关、与 $\theta$ 无关（在积分时 $\hat{\boldsymbol\beta}$ 视为常量），于是

$\int \exp\Big\{-\tfrac12(\boldsymbol\beta-\hat{\boldsymbol\beta})'(X'V^{-1}X)(\boldsymbol\beta-\hat{\boldsymbol\beta})\Big\}\,d\boldsymbol\beta = (2\pi)^{p/2}\,|X'V^{-1}X|^{-1/2}$

整理得到 REML 对数似然：

23.4 Profiling：先消掉 $\boldsymbol\beta$ 与残差尺度

输出里的 Residual Variance Method: Profile 指的是：SAS 把残差方差"profile"出去——即把 $V$ 写成 $\sigma^2 \tilde V(\tilde\theta)$，对给定的 $\tilde\theta$ 解析地求出 $\hat\sigma^2$，再代回目标函数，使优化只在降维后的参数空间进行。这既减少迭代维数，也改善数值稳定性。

$\hat{\boldsymbol\beta}$ 的处理同理：它总是取 GLS 闭式解，因此目标函数实际上只依赖 $\theta$。

23.5 Score 与 Fisher 信息的推导

记

$C = (X'V^{-1}X)^{-1},\qquad P = V^{-1} - V^{-1}XCX'V^{-1},\qquad V_h = \dfrac{\partial V}{\partial\theta_h}$

对 $-2\ell_R$ 逐项求偏导（利用 $\partial V^{-1}/\partial\theta_h = -V^{-1}V_hV^{-1}$）：

关键恒等式：$V^{-1}r = V^{-1}(y-X\hat{\boldsymbol\beta}) = \big(V^{-1} - V^{-1}XCX'V^{-1}\big)y = Py$，因此

$r'V^{-1}V_hV^{-1}r = y'PV_hPy$

又 $\operatorname{tr}(PV_h) = \operatorname{tr}(V^{-1}V_h) - \operatorname{tr}\!\big(CX'V^{-1}V_hV^{-1}X\big)$，三项合并即得：

实现校验：我把解析的 $\mathcal I$ 与用有限差分算出的 $\tfrac12\nabla^2(-2\ell_R)$ 逐一比对，最大相对偏差约 4%（有限差分在目标函数量级为 548、步长 $10^{-4}$ 时的正常误差）。两者一致说明上面的推导与代码实现是自洽的。

23.6 迭代与收敛判据

Fisher scoring / Newton-Raphson 的更新格式为

$\theta^{(t+1)} = \theta^{(t)} + \mathcal{I}\big(\theta^{(t)}\big)^{-1}\, s\big(\theta^{(t)}\big)$

实际实现还要加线搜索：若整步不能使目标函数下降，就把步长折半——真实输出的 Iteration 1 显示 Evaluations = 2，正是这一步多算了一次似然。

23.7 以 CS 为例把导数落到实处

CS 结构下 $\Sigma = \sigma_{cs}J + \sigma^2_e I$，于是

$\dfrac{\partial\Sigma}{\partial\sigma_{cs}} = J,\qquad \dfrac{\partial\Sigma}{\partial\sigma^2_e} = I$

对受试者 $i$ 的分块 $V_i = S_i'\Sigma S_i$（$S_i$ 为选择矩阵），有 $\partial V_i/\partial\theta_h = S_i'(\partial\Sigma/\partial\theta_h)S_i$。代回 23.5 的两条公式即可——本例 $\mathcal I$ 为 $2\times2$：

I = [ 7.84e+00 3.34e+00 ; 3.34e+00 1.997e+01 ]  （自研引擎输出，log 尺度）

求逆并经 delta 法换算回原尺度后得到 SE(CS) = 28.33、SE(Residual) = 14.03，与真实输出置信区间反推的 ≈25、≈13.8 量级一致。

23.8 拟合统计量的 n 口径（数值演算）

用真实输出的数值把三个指标完整算一遍（$d=2$，$N=80$，$\mathrm{rank}(X)=11$，有效受试者 $=36$，$-2\ell_R = 545.71132552$）：

三个公式、三个不同的 $n$（—、69、36）。SAS 文档对此有明确规定：AICC 采用 SAS 6 口径（REML 下 $n = N - \mathrm{rank}(X)$），BIC 采用 Dimensions 表中的有效受试者数。跨软件复核时若统一用一个 $n$，AIC/BIC 必然对不上。

23.9 Kenward-Roger：为什么要调，以及怎么调

动机。$\hat{\boldsymbol\beta}$ 虽是无偏的，但 $\hat\Phi = (X'\hat V^{-1}X)^{-1}$ 作为 $\mathrm{Var}(\hat{\boldsymbol\beta})$ 的估计是有偏的——它把 $\hat\theta$ 当作已知的 $\theta$，忽略了估计协方差参数带来的额外变异。小样本下这会低估标准误、使推断过于乐观。Kackar–Harville 与 Kenward–Roger 的思路是对 $\hat\Phi$ 做泰勒展开修正。

沿用记号 $\Phi = (X'V^{-1}X)^{-1}$，并定义（对协方差参数 $\theta_h$）：

$P_h = X'\dfrac{\partial V^{-1}}{\partial\theta_h}X,\qquad Q_{hj} = X'\dfrac{\partial V^{-1}}{\partial\theta_h}V\dfrac{\partial V^{-1}}{\partial\theta_j}X,\qquad R_{hj} = X'V^{-1}\dfrac{\partial^2 V}{\partial\theta_h\partial\theta_j}V^{-1}X$

令 $W = \mathcal I^{-1}$（$\hat\theta$ 的渐近协方差矩阵）。Kenward-Roger 调整后的协方差矩阵为

分母自由度。设检验 $H_0: K\boldsymbol\beta = 0$，$K$ 为 $c\times p$ 对比矩阵。令 $M = K'(K\Phi K)^{-1}K$，则

$A_1 = \sum_{h,j} W_{hj}\,\mathrm{tr}(M\Phi P_h\Phi)\,\mathrm{tr}(M\Phi P_j\Phi),\qquad A_2 = \sum_{h,j} W_{hj}\,\mathrm{tr}(M\Phi P_h\Phi M\Phi P_j\Phi)$

$B = \dfrac{A_1 + 6A_2}{2c},\qquad g = \dfrac{(c+1)A_1 - (c+4)A_2}{(c+2)A_2},\qquad c_i = \dfrac{g_i}{3c + 2(1-g)}$

$E^* = \Big(1-\dfrac{A_2}{c}\Big)^{-1},\qquad V^* = \dfrac{2}{c}\cdot\dfrac{1 + c_1 B}{(1-c_2 B)^2(1-c_3 B)},\qquad \rho = \dfrac{V^*}{2(E^*)^2}$

其中 $(c_1,c_2,c_3) = \big(g,\ c-g,\ c+2-g\big)\big/\big(3c+2(1-g)\big)$。标量对比时 $c=1$，$F = t^2$，SAS 报告 $t = \hat\psi/\mathrm{SE}_{KR}$、自由度取 $m$。

24.1 先建立因果顺序：一张决策树

这五条提示不是并列的。有些是"模型根本没拟合完"，有些是"拟合完了但结果不可信"，还有些只是需要记录在案的提示。分清层级，才能决定下一步动作。

24.2 Convergence criteria met but final Hessian is not positive definite

后果。协方差参数的标准误来自 $\mathcal I^{-1}$；Hessian 非正定意味着这个矩阵不可信，进而不支持可靠的标准误、CI 和 KR 自由度。

修改策略（按优先级）：

24.3 Asymptotic variance matrix of covariance parameter estimates has been found to be singular and a generalized inverse was used

它在说什么。$\mathcal I$（或其估计）不可逆，SAS 改用广义逆（Moore–Penrose 伪逆）继续计算。这意味着协方差参数之间存在线性依赖——模型里有"多余"的协方差参数。

后果。协方差参数的标准误不可靠。更要紧的是：Kenward-Roger 的计算依赖 $W=\mathcal I^{-1}$，广义逆会改变 $W$，从而影响 $\hat\Phi_A$ 与自由度 $m$。所以这条提示不是"只影响协方差参数那张表"，它会一路传导到主终点的 SE、CI 和 P 值。

修改策略：与 24.2 相同的降复杂度路径；去掉不必要的 GROUP=；检查是否有冗余的协方差结构（例如同时指定 RANDOM 与 REPEATED 造成重复建模）。

24.4 Covariance parameters with zero variance do not contribute to degrees of freedom computed by DDFM=KENWARDROGER

它在说什么。有若干个协方差参数被估成 0（落在边界上）。Kenward-Roger 的自由度公式要对所有协方差参数求和，但零方差参数不携带不确定性信息，因此不计入自由度计算。

如何确认。打开 CovParms，看是否有 Estimate 恰为 0 或极其接近 0 的行。零方差常见于：随机效应方差被估为 0（说明该层变异不存在）、或某访视方差为 0。

处理建议：

• 记录：在 QC 日志与 CSR 脚注中写明"有 $k$ 个零方差协方差参数未计入 KR 自由度"。
• 评估：零方差是否意味着结构选得过大？若是，可考虑更简洁的结构（但这属于统计师决策）。
• 不要为了让这条提示消失而随意加约束或改结构——改变自由度计算基础会影响推断结论。

24.5 Stopped because of too many likelihood evaluations

它在说什么。似然函数的评估次数达到 MAXFUNC=（默认 150）。注意：这与"迭代次数"不是一回事——每次迭代内部若做线搜索（步长折半），会消耗多次评估。

本指南的实跑案例正是这一条。自研引擎在相同骨架数据上拟合 UN（15 参数）时：

UN 迭代史（未收敛现场）：
    iter  0   -2 Res Log Like = 584.77064213
    iter  1   -2 Res Log Like = 584.77064213
    iter  2   -2 Res Log Like = 578.88065859
    iter  3   -2 Res Log Like = 572.17005138
    ...
    iter 12   -2 Res Log Like = 540.99523007
    iter 13   -2 Res Log Like = 540.81729794
    → too_many_likelihood_evaluations

注意这个形态：目标函数仍在持续下降（584.77 → 540.82），但下降速度越来越慢。这正是"似然面平坦 + 参数过多"的典型表现——不是算法坏了，而是数据提供的信息不足以支撑 15 个协方差参数。对比 CS 结构 5 次迭代就收敛，差距一目了然。

修改策略：

24.6 Did not converge

它在说什么。迭代次数达到 MAXITER=（默认 50）而收敛判据仍未满足。与 24.5 的区别在于"卡在哪个计数器上"。

典型形态：目标函数在若干次迭代后来回震荡或下降极其缓慢。常见原因是似然面存在狭长的"山谷"（参数间高度相关），或者参数落在边界附近反复挣扎。

修改策略：与 24.5 基本一致，但额外建议：

• 先用 parms 语句从多个初值试探，判断是"初值问题"还是"结构问题"。
• 若多个合理初值都失败，基本可以确定是数据不支持该结构，应直接走回退链，而不是继续堆 MAXITER。
• 检查是否有近乎完全共线的固定效应（例如某访视只有一个治疗组有数据）。

24.7 贯穿案例：从 UN 失败到 CS 的生产级回退

把上面几条串起来，看真实项目里是怎么处理的（代码已掩蔽：

/* 步骤 1：跑之前先清空，防止不收敛时残留旧结果 */
proc datasets lib=work nolist;
    delete lsm diffs;
quit;

/* 步骤 2：先探 UN，只取 ConvergenceStatus */
ods select none;
ods output ConvergenceStatus=conv_un;
proc mixed data=sub_ana(where=(chg>.)) method=reml;
    class usubjid trt01pn avisitn;
    model chg = base trt01pn avisitn trt01pn*avisitn / ddfm=kr residual;
    repeated avisitn / subject=usubjid type=un;
run;
ods select all;

/* 步骤 3：用 reason 精确判定，而不是看有没有 ERROR */
proc sql noprint;
    select count(*) into :fail_cnt from conv_un
    where reason ne "Convergence criteria met.";
quit;

/* 步骤 4：失败则遍历回退链，记录各结构 AIC，取最小者重跑 */
%let methlist=TOEPH|ARH(1)|CS|CSH|SP(POW)(avisitn_n);

24.8 交付前的最小检查清单