统计分析：从生存曲线到 MMRM

你已经会读 SAP、会审 TLF、知道 HR 与 p 值意味着什么。这一章只做一件事：把 PROC SGPLOT / LIFETEST / PHREG / MIXED / FREQ 换成 R 的等价物，并且让你逐行看懂输出——本章每个代码块都在 R 4.5.0 上实跑过，输出原样贴出。

3.1 ggplot2 的图层语法：把 SGPLOT 的心智拆开重装

SAS 里画图是"一句 PROC SGPLOT + 一堆选项"：你告诉 SAS 画什么图（VBAR）、按什么分组（GROUP=）、怎么排版（ODS）。ggplot2 换了个思路——它不是"画图命令"，而是图形语法（grammar of graphics）：一张图由六个可独立替换的层叠起来。

先用一份自造的迷你 ADSL + ADAE 画 AE 系统器官分类（SOC）的分组条形图。数据是模拟的，结构是真的：USUBJID / TRT01P 在 ADSL，AEDECOD / AEBODSYS 在 ADAE，按 USUBJID 合并——和你在 SAS 里做的完全一样。

library(ggplot2)
set.seed(2024)
n <- 120
## 迷你 ADSL：120 例受试者，随机分到三组
adsl <- data.frame(
  USUBJID = sprintf("S%03d", seq_len(n)),
  TRT01P  = factor(rep(c("Placebo", "Drug 20mg", "Drug 50mg"), length.out = n),
                   levels = c("Placebo", "Drug 20mg", "Drug 50mg"))
)
## 迷你 ADAE：214 条 AE 记录，PT 归属 5 个 SOC
soc  <- c("Nervous system", "Gastrointestinal", "Skin", "General", "Respiratory")
pt   <- list(c("Headache", "Dizziness"), c("Nausea", "Diarrhoea"),
             c("Rash", "Pruritus"), c("Fatigue", "Pyrexia"), c("Cough"))
aeterm <- unlist(pt)
aesoc  <- rep(soc, sapply(pt, length))
ae <- data.frame(
  USUBJID  = sample(adsl$USUBJID, 214, replace = TRUE),
  AEDECOD  = sample(aeterm, 214, replace = TRUE)
)
ae$AEBODSYS <- aesoc[match(ae$AEDECOD, aeterm)]
ae <- merge(ae, adsl, by = "USUBJID")   # 等价于 SAS: merge ae adsl; by USUBJID;

## 四层就够出图：data -> mapping -> geom -> theme
p <- ggplot(ae, aes(x = AEBODSYS, fill = TRT01P)) +
  geom_bar(position = "dodge") +
  labs(x = "System Organ Class", y = "Number of AE events",
       title = "AE count by SOC and treatment") +
  theme_minimal(base_size = 12) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggsave("fig31_ae_soc.png", p, width = 7, height = 4.2, dpi = 150)
cat("png written:", file.exists("fig31_ae_soc.png"), "\n")
print(table(ae$AEBODSYS, ae$TRT01P))
cat("subjects with >=1 AE:", length(unique(ae$USUBJID)), "of", n, "\n")

# 实跑捕获（R 4.5.0，已去空行）
png written: TRUE
                   Placebo Drug 20mg Drug 50mg
  Gastrointestinal      17        19        16
  General               19        14        20
  Nervous system        15        19        17
  Respiratory           11         6         3
  Skin                  12        12        14
subjects with >=1 AE: 102 of 120

再叠两层看看威力——scale 自定义颜色与坐标轴余量，facet 按治疗组拆面板。同样的数据，换两层就换了一种表达：

## scale 层：自定义填充色 + 去掉 y 轴下方留白；facet 层：按治疗组拆面板
adsl$SEX <- factor(sample(c("F", "M"), n, replace = TRUE))
ae <- merge(ae, adsl[, c("USUBJID", "SEX")], by = "USUBJID")
p2 <- ggplot(ae, aes(x = AEBODSYS, fill = SEX)) +
  geom_bar(position = "dodge") +
  facet_wrap(~ TRT01P) +
  scale_fill_manual(values = c(F = "#8172B3", M = "#64B5CD"), name = "Sex") +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
  labs(x = "System Organ Class", y = "Number of AE events") +
  theme_minimal(base_size = 11) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggsave("fig31b_ae_facet.png", p2, width = 9, height = 3.6, dpi = 150)
cat("facet png written:", file.exists("fig31b_ae_facet.png"), "\n")
print(table(ae$TRT01P, ae$SEX))

# 实跑捕获（R 4.5.0，已去空行）
facet png written: TRUE
             F  M
  Placebo   32 42
  Drug 20mg 38 32
  Drug 50mg 31 39

3.2 Surv 对象：生存数据在 R 里的最小结构

SAS 的 PROC LIFETEST 吃两个变量：TIME= 和 STATUS=。R 的 survival 包把这两列打包成一个特殊矩阵——Surv 对象，删失的行在打印时带 + 号。

library(survival)
lung <- survival::lung
## Surv(时间, 事件指示) -> 一个 228 x 2 的特殊矩阵
s <- Surv(lung$time, lung$status == 2)
cat("class:", class(s), "| type:", attr(s, "type"),
    "| dim:", paste(dim(s), collapse = " x "), "\n")
print(head(s, 10))
cat("--- event flag column, first 10 ---\n")
print(head(s[, "status"], 10))
cat("--- original status coding in lung ---\n")
print(table(lung$status))
cat("--- recoded event counts (1 = death, 0 = censored) ---\n")
print(table(s[, "status"]))

# 实跑捕获（R 4.5.0，已去空行）
class: Surv | type: right | dim: 228 x 2
 [1]  306   455  1010+  210   883  1022+  310   361   218   166
--- event flag column, first 10 ---
 [1] 1 1 0 1 1 0 1 1 1 1
--- original status coding in lung ---
  1   2
 63 165
--- recoded event counts (1 = death, 0 = censored) ---
  0   1
 63 165

3.3 KM 曲线与 log-rank 检验

三件套：survfit() 出 KM 估计与中位生存时间、survdiff() 出 log-rank 检验、survminer::ggsurvplot() 出可直接进报告的图。数据仍用 lung，按性别分组。

library(survival)
lung <- survival::lung
lung$event <- as.integer(lung$status == 2)   # 1 = 死亡，0 = 删失
fit <- survfit(Surv(time, event) ~ sex, data = lung)
print(fit)                                   # n / events / median / 95% CI
cat("--- survival rate at 365 days ---\n")
print(round(summary(fit, times = 365)$surv, 3))
cat("--- log-rank test ---\n")
print(survdiff(Surv(time, event) ~ sex, data = lung))
cat("--- survminer plot ---\n")
library(survminer)
p <- ggsurvplot(fit, data = lung, pval = TRUE, conf.int = TRUE,
                risk.table = TRUE, xlab = "Days", ylab = "Survival probability",
                legend.labs = c("Male", "Female"),
                palette = c("#4C72B0", "#DD8452"), ggtheme = theme_minimal())
ggsave("fig33_km_sex.png", plot = p$plot, width = 7, height = 4.6, dpi = 150)
cat("png written:", file.exists("fig33_km_sex.png"), "\n")

# 实跑捕获（R 4.5.0，已去空行）
Call: survfit(formula = Surv(time, event) ~ sex, data = lung)
        n events median 0.95LCL 0.95UCL
sex=1 138    112    270     212     310
sex=2  90     53    426     348     550
--- survival rate at 365 days ---
[1] 0.336 0.526
--- log-rank test ---
Call:
survdiff(formula = Surv(time, event) ~ sex, data = lung)
        N Observed Expected (O-E)^2/E (O-E)^2/V
sex=1 138      112     91.6      4.55      10.3
sex=2  90       53     73.4      5.68      10.3
 Chisq= 10.3  on 1 degrees of freedom, p= 0.001
--- survminer plot ---
png written: TRUE

3.4 Cox 比例风险模型

coxph() 就是 PROC PHREG。模型：Surv(time, event) ~ age + sex + ph.ecog。

library(survival)
lung <- survival::lung
lung$event <- as.integer(lung$status == 2)
d <- lung[, c("time", "event", "age", "sex", "ph.ecog")]
cat("rows total:", nrow(d), "| after dropping NA:", nrow(na.omit(d)), "\n")
d <- na.omit(d)
fit <- coxph(Surv(time, event) ~ age + sex + ph.ecog, data = d)
s <- summary(fit)
cat("ties method:", fit$method, "\n")
print(round(s$coefficients, 4))
cat("--- hazard ratios (exp(coef)) with 95% CI ---\n")
print(round(s$conf.int, 3))
cat("concordance:", round(s$concordance[1], 3),
    "| likelihood ratio p =", signif(s$logtest[3], 3), "\n")
cat("n =", s$n, "| events =", s$nevent, "\n")

# 实跑捕获（R 4.5.0）
rows total: 228 | after dropping NA: 227
ties method: efron
           coef exp(coef) se(coef)       z Pr(>|z|)
age      0.0111    1.0111   0.0093  1.1942   0.2324
sex     -0.5526    0.5754   0.1677 -3.2945   0.0010
ph.ecog  0.4637    1.5900   0.1136  4.0829   0.0000
--- hazard ratios (exp(coef)) with 95% CI ---
        exp(coef) exp(-coef) lower .95 upper .95
age         1.011      0.989     0.993     1.030
sex         0.575      1.738     0.414     0.799
ph.ecog     1.590      0.629     1.273     1.986
concordance: 0.637 | likelihood ratio p = 1.08e-06
n = 227 | events = 164

3.5 MMRM：mmrm 包与 PROC MIXED 的对照

MMRM（混合效应重复测量模型）是纵向疗效分析的主力，R 里的对应物是 mmrm 包——它不是 lme4 的包装，而是基于 TMB 的专用实现，默认就给你 Satterthwaite（或 Kenward-Roger）自由度调整，专为"每个访视点的治疗差"这一估计目标设计。

先造数据：80 例受试者 × 3 个访视（Week 4/8/12），变量 USUBJID / ARM / BASE / AVISIT / AVISITN / FEV1，并在 Week 12 制造 14 条 MAR 式脱落。

set.seed(42)
nsub <- 80
sub <- data.frame(
  USUBJID = sprintf("S%03d", seq_len(nsub)),
  ARM     = factor(rep(c("Placebo", "Active"), length.out = nsub),
                   levels = c("Placebo", "Active")),
  BASE    = round(rnorm(nsub, mean = 50, sd = 10), 1)
)
sub$bi <- round(rnorm(nsub, 0, 4), 2)   # 受试者随机截距：纵向相关性的来源
dat <- data.frame(
  USUBJID = rep(sub$USUBJID, each = 3),
  ARM     = factor(rep(sub$ARM, each = 3), levels = levels(sub$ARM)),
  BASE    = rep(sub$BASE, each = 3),
  bi      = rep(sub$bi, each = 3),
  AVISITN = rep(1:3, times = nsub),
  AVISIT  = factor(rep(c("Week 4", "Week 8", "Week 12"), times = nsub),
                   levels = c("Week 4", "Week 8", "Week 12"))
)
dat$FEV1 <- round(dat$BASE + 3 * dat$AVISITN +
                  ifelse(dat$ARM == "Active", 2.5 * dat$AVISITN, 0) +
                  dat$bi + rnorm(nrow(dat), 0, 3), 1)
## MAR 式脱落：随机删掉 14 条 Week 12 记录
dropout <- sample(which(dat$AVISITN == 3), 14)
dat <- dat[!seq_len(nrow(dat)) %in% dropout,
           c("USUBJID", "ARM", "BASE", "AVISITN", "AVISIT", "FEV1")]
cat("rows:", nrow(dat), "| subjects:", length(unique(dat$USUBJID)), "\n")
print(table(dat$ARM, dat$AVISIT))
cat(paste(capture.output(head(dat, 4)), collapse = "\n"), "\n")

# 实跑捕获（R 4.5.0，已去空行）
rows: 226 | subjects: 80
          Week 4 Week 8 Week 12
  Placebo     40     40      35
  Active      40     40      31
  USUBJID     ARM BASE AVISITN  AVISIT FEV1
1    S001 Placebo 63.7       1  Week 4 72.2
2    S001 Placebo 63.7       2  Week 8 72.5
3    S001 Placebo 63.7       3 Week 12 79.2
4    S002  Active 44.4       1  Week 4 49.8

library(mmrm)
fit <- mmrm(FEV1 ~ ARM + AVISIT + ARM:AVISIT + BASE + us(AVISIT | USUBJID),
            data = dat)
print(fit)

# 实跑捕获（R 4.5.0，已去空行）
mmrm fit
Formula:     FEV1 ~ ARM + AVISIT + ARM:AVISIT + BASE + us(AVISIT | USUBJID)
Data:        dat (used 226 observations from 80 subjects with maximum 3
timepoints)
Covariance:  unstructured (6 variance parameters)
Inference:   REML
Deviance:    1252.894
Coefficients:
            (Intercept)               ARMActive            AVISITWeek 8
              3.6621742               3.0170305               2.8625000
          AVISITWeek 12                    BASE  ARMActive:AVISITWeek 8
              5.6948399               0.9731705               2.9075000
ARMActive:AVISITWeek 12
              5.7788999
Model Inference Optimization:
Converged with code 0 and message: convergence: rel_reduction_of_f <= factr*epsmch

模型跑通只是第一步，你要的三个东西是：协方差矩阵、治疗差的 CI、指定访视点的校正均数差。

cat("--- residual within-subject covariance (US, 3 visits) ---\n")
print(round(VarCorr(fit), 2))
cat("--- treatment-related coefficients with 95% CI ---\n")
tab <- cbind(Estimate = coef(fit), confint(fit))
print(round(tab[grep("ARM", rownames(tab)), , drop = FALSE], 3))
cat("--- adjusted difference at Week 12 via emmeans ---\n")
library(emmeans); emm_options(lmer.df = "satterthwaite")
em <- emmeans(fit, pairwise ~ ARM | AVISIT)
print(subset(em$contrasts, AVISIT == "Week 12"))
cat("--- information criteria ---\n")
cat("logLik:", round(as.numeric(logLik(fit)), 2), "| AIC:", round(AIC(fit), 2),
    "| BIC:", round(BIC(fit), 2), "\n")

# 实跑捕获（R 4.5.0，已去空行）
--- residual within-subject covariance (US, 3 visits) ---
        Week 4 Week 8 Week 12
Week 4   25.06  12.13   16.54
Week 8   12.13  18.71   14.15
Week 12  16.54  14.15   24.77
--- treatment-related coefficients with 95% CI ---
                        Estimate 2.5 % 97.5 %
ARMActive                  3.017 0.788  5.246
ARMActive:AVISITWeek 8     2.908 0.941  4.874
ARMActive:AVISITWeek 12    5.779 3.829  7.729
--- adjusted difference at Week 12 via emmeans ---
 contrast         AVISIT  estimate   SE   df t.ratio p.value
 Placebo - Active Week 12     -8.8 1.16 75.5  -7.551 <0.0001
--- information criteria ---
logLik: -626.45 | AIC: 1264.89 | BIC: 1279.19

3.6 lme4 对照：什么时候用 lmer，什么时候用 mmrm

同一份数据，用 lme4::lmer() 写"随机截距模型"。它更通用、生态更大，但假设更强。

library(lme4)
fit2 <- lmer(FEV1 ~ ARM * AVISIT + BASE + (1 | USUBJID), data = dat)
cat("--- variance components: random intercept => compound symmetry ---\n")
print(VarCorr(fit2))
cat("--- fixed effects ---\n")
print(round(summary(fit2)$coefficients[, 1:2], 3))
cat("REML deviance:", round(deviance(fit2, REML = TRUE), 2),
    "| logLik:", round(as.numeric(logLik(fit2)), 2), "\n")
library(mmrm)
fit1 <- mmrm(FEV1 ~ ARM + AVISIT + ARM:AVISIT + BASE + us(AVISIT | USUBJID), data = dat)
cat("mmrm US logLik:", round(as.numeric(logLik(fit1)), 2),
    "| mmrm US AIC:", round(AIC(fit1), 2), "\n")
cat("lmer RI AIC:", round(AIC(fit2), 2), "(ML vs REML: not directly comparable)\n")
fit2ml <- lmer(FEV1 ~ ARM * AVISIT + BASE + (1 | USUBJID), data = dat, REML = FALSE)
cat("lmer ML AIC:", round(AIC(fit2ml), 2),
    "| mmrm ML AIC:", round(AIC(update(fit1, reml = FALSE)), 2), "\n")

# 实跑捕获（R 4.5.0）
--- variance components: random intercept => compound symmetry ---
 Groups   Name        Std.Dev.
 USUBJID  (Intercept) 3.7517
 Residual             2.9443
--- fixed effects ---
                        Estimate Std. Error
(Intercept)                2.999      2.308
ARMActive                  3.017      1.066
AVISITWeek 8               2.862      0.658
AVISITWeek 12              5.700      0.690
BASE                       0.986      0.043
ARMActive:AVISITWeek 8     2.908      0.931
ARMActive:AVISITWeek 12    5.848      0.998
REML deviance: 1259.42 | logLik: -629.09
mmrm US logLik: -626.45 | mmrm US AIC: 1264.89
lmer RI AIC: 1276.19 (ML vs REML: not directly comparable)
lmer ML AIC: 1277.42 | mmrm ML AIC: 1265.9

3.7 分类数据：卡方与 CMH 分层分析

响应率比较：先粗分析（2×2 卡方），再按中心分层（Mantel-Haenszel）。下面四个中心的合计正好等于粗表，这样你能直接看出分层有没有改变结论。

tab <- matrix(c(42, 18, 30, 35), nrow = 2, byrow = TRUE,
              dimnames = list(Treatment = c("Active", "Placebo"),
                              Response  = c("Yes", "No")))
print(tab)
print(chisq.test(tab))
cat("crude OR:", round((42 * 35) / (18 * 30), 3), "\n")

# 实跑捕获（R 4.5.0，已去空行）
         Response
Treatment Yes No
  Active   42 18
  Placebo  30 35
	Pearson's Chi-squared test with Yates' continuity correction
data:  tab
X-squared = 6.3209, df = 1, p-value = 0.01193
crude OR: 2.722

## 2 x 2 x K 数组：第三维是分层变量（中心）
cen <- array(dim = c(2, 2, 4),
             dimnames = list(Treatment = c("Active", "Placebo"),
                             Response  = c("Yes", "No"),
                             Center    = paste("Site", 1:4)))
vals <- list(c(12, 4, 9, 11), c(9, 3, 8, 10), c(11, 6, 7, 6), c(10, 5, 6, 8))
for (k in 1:4) cen[, , k] <- matrix(vals[[k]], 2, byrow = TRUE)
cat("--- site margins sum back to the crude table ---\n")
print(apply(cen, c(1, 2), sum))
print(mantelhaen.test(cen))

# 实跑捕获（R 4.5.0，已去空行）
--- site margins sum back to the crude table ---
         Response
Treatment Yes No
  Active   42 18
  Placebo  30 35
	Mantel-Haenszel chi-squared test with continuity correction
data:  cen
Mantel-Haenszel X-squared = 6.1387, df = 1, p-value = 0.01323
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
 1.305296 5.774277
sample estimates:
common odds ratio
         2.745385

3.8 森林图：forestploter 画出亚组 HR

森林图的数据框要有一个留白列（这里名字就是一个空格 ` `），图形会画在那一列的位置上——这是 forestploter 的核心机制。

library(forestploter)
library(grid)
library(ggplot2)
df <- data.frame(
  Subgroup = c("Overall", "  Age < 65", "  Age >= 65", "  Male", "  Female"),
  N        = c(420, 250, 170, 240, 180),
  `HR (95% CI)` = c("0.72 (0.55-0.94)", "0.65 (0.46-0.92)", "0.84 (0.55-1.28)",
                    "0.70 (0.50-0.98)", "0.76 (0.48-1.20)"),
  ` `      = rep("", 5),
  check.names = FALSE, stringsAsFactors = FALSE)
est  <- c(0.72, 0.65, 0.84, 0.70, 0.76)
low  <- c(0.55, 0.46, 0.55, 0.50, 0.48)
high <- c(0.94, 0.92, 1.28, 0.98, 1.20)
tm <- forest_theme(base_size = 10, refline_gp = gpar(col = "grey40"),
                   ci_Theight = 0.01, legend.position = "bottom")
p <- forest(df, est = est, lower = low, upper = high,
            ci_column = 4, ref_line = 1, xlim = c(0.4, 1.4), theme = tm)
ggsave("fig38_forest.png", p, width = 8, height = 3.7, dpi = 150)
cat("class(p):", paste(class(p), collapse = "/"), "| png written:",
    file.exists("fig38_forest.png"), "\n")
print(df[, c("Subgroup", "N", "HR (95% CI)")])

# 实跑捕获（R 4.5.0）
class(p): forestplot/gtable/gTree/grob/gDesc | png written: TRUE
     Subgroup   N      HR (95% CI)
1     Overall 420 0.72 (0.55-0.94)
2    Age < 65 250 0.65 (0.46-0.92)
3   Age >= 65 170 0.84 (0.55-1.28)
4        Male 240 0.70 (0.50-0.98)
5      Female 180 0.76 (0.48-1.20)

3.9 期中分析与成组序贯边界：gsDesign

期中分析的数学核心是 α 消耗（alpha spending）：你多看一次数据，就必须为这次"看"预付一部分 I 类错误预算。gsDesign 直接算出边界。

library(gsDesign)
## k = 2 次分析（1 次期中 + 1 次最终），test.type = 1 为单侧优效设计
## alpha/beta/timing 均为默认值，这里显式写出便于对照 SAP
d <- gsDesign(k = 2, test.type = 1, alpha = 0.025, beta = 0.1, timing = 1)
print(d)

# 实跑捕获（R 4.5.0，已去空行）
One-sided group sequential design with
90 % power and 2.5 % Type I Error.
           Sample
            Size
  Analysis Ratio*  Z   Nominal p Spend
         1  0.504 2.75    0.0030 0.003
         2  1.009 1.98    0.0238 0.022
     Total                       0.0250
++ alpha spending:
 Hwang-Shih-DeCani spending function with gamma = -4.
* Sample size ratio compared to fixed design with no interim
Boundary crossing probabilities and expected sample size
assume any cross stops the trial
Upper boundary (power or Type I Error)
          Analysis
   Theta      1      2 Total   E{N}
  0.0000 0.0030 0.0220 0.025 1.0072
  3.2415 0.3271 0.5729 0.900 0.8437

cat("--- bounds on the Z scale ---\n")
print(round(d$upper$bound, 4))
cat("--- information fraction / nominal p ---\n")
print(round(cbind(timing = d$timing, z = d$upper$bound,
                  p_nominal = pnorm(d$upper$bound, lower.tail = FALSE)), 4))
cat("information at analysis 1 / 2:", round(d$n.I, 3), "\n")

# 实跑捕获（R 4.5.0）
--- bounds on the Z scale ---
[1] 2.7500 1.9811
--- information fraction / nominal p ---
     timing      z p_nominal
[1,]    0.5 2.7500    0.0030
[2,]    1.0 1.9811    0.0238
information at analysis 1 / 2: 0.504 1.009

章末资源

本章练习

1. （易）用 lung 把 KM 按 ECOG 体能状态分层：造 event 列，取 ph.ecog %in% c(0, 1) 的子集（ph.ecog = 3 只有 1 例，画出来没有意义），做 factor 后 survfit() + survdiff()。参考结果（已实跑）：ECOG 0 组 n=63、中位 394 天（348–574）；ECOG 1 组 n=113、中位 306 天（268–429）；log-rank Chisq = 3.5，p = 0.06。
2. （中）在 3.4 的 Cox 模型里加交互项 sex * ph.ecog，用 anova(fit0, fit1) 做似然比检验，并解释"加了交互项后主效应系数的含义变了"。参考结果（已实跑）：交互项 exp(coef) = 1.178、p = 0.51；LR 检验 Chisq = 0.436、Df = 1、p = 0.509 ⇒ 无证据支持交互。注意此时 sex 的系数（−0.730）已不再是"总体性别效应"，而是 ph.ecog = 0 时的性别效应。
3. （难）把 3.5 的 MMRM 换成不同协方差结构，用 ML（reml = FALSE）比较 AIC/BIC，并检查 Week 12 的校正差是否稳健：us() / cs() / ar1() / toeph()。参考结果（已实跑）：AIC 分别 1265.90 / 1263.42 / 1280.07 / 1265.71，BIC 1280.19 / 1268.19 / 1284.83 / 1277.62 ⇒ CS 在本数据上最优（参数少、拟合几乎不差）；Week 12 校正差在 8.80–8.96 之间，结论对协方差结构不敏感。想一想：为什么 REML 下不能这样比 AIC？（提示：REML 似然是"固定效应被投影掉"之后的似然，只在固定效应完全相同时可比。）

下一章：统计结果怎么变成能递交的表格与图——rtables / flextable / Quarto 搭一条 TLF 流水线，把本章的 KM 表、Cox 表、MMRM 表一次性产出并双编程比对（第 4 章）。