The V Lab
适用对象医学、公共卫生、心理学与健康服务研究学习者
学习时长约 180–240 分钟
先修要求线性回归、置信区间、矩阵与基本 R 语法

贯穿案例 360 名高血压管理项目参与者随机接受常规照护或干预，在第 0、3、6、9、12 月测量收缩压（SBP）。每人的起点和变化速度不同，同一个人的相邻测量还存在 AR(1) 残差相关。随访概率只依赖已经观测到的上一访视 SBP 和随机组，因此在生成机制下属于“给定已观测历史”的 MAR。所有记录均为固定种子的模拟资料，不是临床证据。

纵向模型不是因果识别器 模型可以正确表达同一个人的相关测量，却不会自动消除未测混杂、选择性入组、干预污染或 MNAR 失访。随机分配支持本案例的组间因果比较；对观察性暴露，仍须单独陈述识别假设。

如何使用本教程

推荐沿着下面的研究链条学习：

研究问题与 estimand → 数据结构和质量 → 轨迹可视化 → 时间函数与相关结构 → 边际模型／混合模型 → 模型比较和诊断 → 缺失与敏感性分析 → 报告

学习目标

完成本教程后，你应能够：

  • 区分纵向、重复横断面和一般聚类资料；
  • 在 long/wide 格式之间转换，并审计重复访视、时间顺序和失访模式；
  • 选择连续、分类、二次或样条时间表达；
  • 解释独立、复合对称（CS）、AR(1)、非结构化协方差与随机效应；
  • 说明普通 OLS 为何会给出错误的不确定性；
  • 区分 population-averaged 与 subject-specific estimand；
  • 用 gls()、MMRM 和 lme() 分析连续纵向结局；
  • 从模型矩阵标准化得到各组访视均值、12 月组差和变化差；
  • 区分总体轨迹预测与已有参与者的条件预测；
  • 识别 MCAR、MAR、MNAR，说明似然分析、MI、IPW 和联合模型的边界；
  • 处理 time-varying covariate 的 within/between 分解；
  • 为纵向二元或计数结局选择 GEE 或 GLMM，并透明报告限制。

1 研究设计与数据结构

1.1 纵向不等于“有时间变量”

结构 同一人是否重复 典型目标 主要依赖来源
纵向队列／试验 是 个体内变化、组间轨迹差 人内相关、失访
重复横断面 否，波次抽不同人 总体均值随日历时间变化 波次内抽样／聚类
单次聚类资料 否 医院、学校或社区差异 群组内相关
密集纵向资料 是，频率很高 短时动态、滞后关系 自相关、昼夜周期、不规则间隔

“第 0、3、6 月各调查 500 人”若每次人员不同，是重复横断面；不能计算个人变化。相反，每人只有两次测量也已是纵向资料。分析单位不是数据行，而是研究问题、随机化单位和相关结构共同决定的。

1.2 Long 格式审计

audit_summary <- data.frame(
  项目 = c("参与者数", "计划记录数", "实际观测数", "重复 id-time 键", "时间范围（月）"),
  结果 = c(
    length(unique(long_data$id_num)), nrow(long_data), sum(long_data$observed),
    sum(duplicated(long_data[c("id_num", "time")])),
    paste(range(long_data$time), collapse = "–")
  )
)
knitr::kable(audit_summary, caption = "纵向 long 数据的最小质量审计")
纵向 long 数据的最小质量审计
项目 结果
参与者数 360
计划记录数 1800
实际观测数 1642
重复 id-time 键 0
时间范围（月） 0–12

应进一步核对：基线是否早于干预、时间单位是否一致、同一次访视是否出现两条记录、结局单位是否改变、访视窗口如何映射到计划时间、缺失码是否误作数值、以及参与者是否在中途改变随机组。若实际时间偏离较大，可同时保存 scheduled_time 与 actual_time。

1.3 Wide 格式与往返转换

wide_sbp <- reshape(
  long_data[c("id_num", "time", "sbp")], idvar = "id_num",
  timevar = "time", direction = "wide"
)
names(wide_sbp) <- sub("sbp\\.", "SBP_月", names(wide_sbp))
knitr::kable(head(wide_sbp, 8), digits = 1, caption = "前 8 人的 wide 格式 SBP")
前 8 人的 wide 格式 SBP
id_num SBP_月0 SBP_月3 SBP_月6 SBP_月9 SBP_月12
1 1 146 143 141 138 135
6 2 151 149 147 143 150
11 3 151 155 149 148 156
16 4 129 128 135 128 128
21 5 140 141 136 131 131
26 6 139 140 129 135 135
31 7 142 147 151 142 140
36 8 142 143 138 136 135 
round_trip <- reshape(
  wide_sbp, idvar = "id_num", varying = names(wide_sbp)[-1],
  v.names = "sbp", timevar = "time_label", direction = "long"
)
stopifnot(nrow(round_trip) == n * length(times))

wide 格式适合人工浏览、制作基线至终点变化量；long 格式通常更适合混合模型和不等次数随访。转换前必须验证 id × time 唯一，否则 reshape() 可能隐藏重复记录。

1.4 访视保留与单调模式

retention <- aggregate(observed ~ time, long_data, mean)
names(retention) <- c("月", "保留比例")
retention$保留百分比 <- sprintf("%.1f%%", 100 * retention$保留比例)

pattern_string <- apply(observed_matrix, 1, paste0, collapse = "")
pattern_table <- sort(table(pattern_string), decreasing = TRUE)
pattern_display <- data.frame(
  模式 = names(pattern_table), 人数 = as.integer(pattern_table),
  百分比 = sprintf("%.1f%%", 100 * as.integer(pattern_table) / n)
)
knitr::kable(retention[c("月", "保留百分比")], caption = "各计划访视的资料保留率")
各计划访视的资料保留率
月 保留百分比
0 100.0%
3 94.4%
6 90.8%
9 87.5%
12 83.3% 
knitr::kable(pattern_display, caption = "1=观测、0=缺失的单调访视模式")
1=观测、0=缺失的单调访视模式
模式 人数 百分比
TRUETRUETRUETRUETRUE 300 83.3%
TRUEFALSEFALSEFALSEFALSE 20 5.6%
TRUETRUETRUETRUEFALSE 15 4.2%
TRUETRUEFALSEFALSEFALSE 13 3.6%
TRUETRUETRUEFALSEFALSE 12 3.3%

本例一旦退出，之后都缺失；真实资料也可能出现间歇缺失。不要只报告“总缺失百分比”：各组、各访视、退出前结局和退出原因都与可接受的缺失假设有关。

2 先看轨迹，再拟合模型

2.1 个体 spaghetti plot

set_long_plot_font()
set.seed(20261006)
shown_ids <- sort(sample(id, 40))
shown <- analysis_data[analysis_data$id_num %in% shown_ids, ]
op <- par(mar = c(4.4, 4.4, 1.0, 1.0))
plot(NA, xlim = range(times), ylim = range(shown$sbp), xlab = "月份", ylab = "SBP（mmHg）")
for (one_id in shown_ids) {
  d <- shown[shown$id_num == one_id, ]
  col <- if (d$program_num[1] == 1) palette_long[["orange"]] else palette_long[["blue"]]
  lines(d$time, d$sbp, type = "o", pch = 16, cex = 0.45,
        col = grDevices::adjustcolor(col, 0.38))
}
legend("top", inset = -0.02, horiz = TRUE, bty = "n",
       legend = c("常规照护", "干预"), col = c(palette_long[["blue"]], palette_long[["orange"]]),
       lty = 1, pch = 16)
40名参与者的收缩压随月份变化的折线图，颜色区分常规照护和干预。

抽取 40 名参与者的观测 SBP 轨迹；断线表示退出后没有资料。 

par(op)

图中既有组间趋势，也有明显的个体起点和斜率差异。Spaghetti plot 的目的不是靠肉眼“证明显著”，而是发现不合理跳变、天花板/地板、非线性、异常访视和可能需要的随机斜率。

2.2 各组均值与不确定性

mean_by_visit <- do.call(rbind, lapply(split(analysis_data, list(analysis_data$program, analysis_data$time)), function(d) {
  n_obs <- nrow(d)
  data.frame(program = d$program[1], time = d$time[1], n = n_obs,
             mean = mean(d$sbp), se = sd(d$sbp) / sqrt(n_obs))
}))
mean_by_visit <- mean_by_visit[order(mean_by_visit$program, mean_by_visit$time), ]
mean_by_visit$lower <- with(mean_by_visit, mean - qt(0.975, n - 1) * se)
mean_by_visit$upper <- with(mean_by_visit, mean + qt(0.975, n - 1) * se)
set_long_plot_font()
op <- par(mar = c(4.4, 4.4, 2.8, 1.0), xpd = NA)
plot(NA, xlim = range(times), ylim = range(c(mean_by_visit$lower, mean_by_visit$upper)) + c(-1, 2),
     xlab = "月份", ylab = "观测 SBP 均值（mmHg）")
for (g in levels(long_data$program)) {
  d <- mean_by_visit[mean_by_visit$program == g, ]
  col <- if (g == "Intervention") palette_long[["orange"]] else palette_long[["blue"]]
  lines(d$time, d$mean, type = "o", pch = 16, lwd = 2, col = col)
  arrows(d$time, d$lower, d$time, d$upper, angle = 90, code = 3, length = 0.04, col = col)
}
legend("top", inset = -0.13, horiz = TRUE, bty = "n",
       legend = c("常规照护", "干预"), col = c(palette_long[["blue"]], palette_long[["orange"]]),
       lty = 1, lwd = 2, pch = 16)
常规照护和干预两组的收缩压均值随月份变化，带置信区间误差棒。

按随机组和访视计算的观测 SBP 均值及 95% t 置信区间。 

par(op)

这些是各访视仍被观测者的原始均值，不一定代表最初随机入组总体。误差棒也没有校正重复测量相关性；它们是描述，不是主要纵向检验。

2.3 基线到 12 月变化

complete_12 <- wide_sbp[!is.na(wide_sbp$SBP_月0) & !is.na(wide_sbp$SBP_月12), ]
complete_12$change <- complete_12$SBP_月12 - complete_12$SBP_月0
complete_12$program <- factor(program_num[complete_12$id_num], levels = 0:1,
                              labels = c("常规照护", "干预"))
set_long_plot_font()
op <- par(mar = c(4.4, 4.4, 1.0, 1.0))
boxplot(change ~ program, data = complete_12, col = c("#DCEAF7", "#FBE4BE"),
        ylab = "12 月 - 基线 SBP（mmHg）", xlab = "随机组", outline = FALSE)
stripchart(change ~ program, data = complete_12, vertical = TRUE, method = "jitter",
           pch = 16, cex = 0.45,
           col = grDevices::adjustcolor(palette_long[["navy"]], 0.35), add = TRUE)
abline(h = 0, lty = 2, col = palette_long[["gray"]])
两组基线到12月收缩压变化量的箱线图和抖动散点。

有基线和 12 月资料者的个体变化量；负值表示 SBP 下降。 

par(op)

只分析完整 12 月资料会丢掉中间轨迹，并在退出与结局有关时改变目标总体。变化量也会把基线测量误差带入响应。它可作为描述或预设敏感性分析，但不是纵向模型的默认替代品。

3 时间与协方差怎么编码

3.1 连续、分类、二次与样条时间

时间表达 模型含义 优点 风险／适用场景
连续线性 time 每月恒定变化 简洁、功效高 错过弯曲或平台期
分类 factor(time) 每个访视独立均值 少依赖形状假设 参数多，不能自然内插
二次 time + I(time^2) 平滑弯曲 本例生成机制匹配 边界外外推危险
样条 splines::ns(time, df=3) 分段平滑 灵活 自由度和 knots 应预设

若模型含 program * time，programIntervention 是 time=0 的组差；programIntervention:time 是两组每月线性斜率之差。把时间中心化到 12 月会让组主项变为 12 月差，但不会改变拟合值。二次模型中任一时点的瞬时斜率还包含 2βt2t2\beta_{t^2}t。

基线处理取决于设计和 estimand：可将基线作为纵向响应的一部分并约束随机试验的共同基线，也可分析随访结局并调整个体基线（cLDA/ANCOVA 思路）。不能既把同一基线作为响应又机械地复制为每行协变量，却不说明模型约束。

3.2 相关结构的语言

对同一个体 ii 的响应向量 𝐘i\mathbf Y_i，常见候选包括：

  • 独立：条件协方差矩阵为对角；纵向资料通常不可信；
  • CS／交换型：任意两次相关相同，适合无明显时间衰减的情形；
  • AR(1)：相邻访视相关为 ρ\rho，间隔 kk 个计划访视为 ρk\rho^k；等间隔假设很重要；
  • 非结构化：每个方差和协方差单独估计，灵活但耗参数；
  • 随机截距／斜率：通过个体轨迹异质性诱导相关，可再叠加残差 AR(1)。
complete_all <- wide_sbp[complete.cases(wide_sbp), -1]
emp_cor <- cor(complete_all)
set_long_plot_font()
op <- par(mar = c(5.3, 5.3, 1.0, 1.0))
image(seq_along(times), seq_along(times), emp_cor[nrow(emp_cor):1, ],
      axes = FALSE, xlab = "访视月", ylab = "访视月",
      col = colorRampPalette(c("white", palette_long[["sky"]], palette_long[["navy"]]))(40), zlim = c(0, 1))
axis(1, at = seq_along(times), labels = times, las = 1, padj = 0.7)
axis(2, at = seq_along(times), labels = rev(times), las = 1)
for (i in seq_along(times)) for (j in seq_along(times)) {
  text(j, length(times) - i + 1, sprintf("%.2f", emp_cor[i, j]), cex = 0.8)
}
五次访视收缩压经验相关系数的蓝色热图。

完成全部五次访视者的经验相关矩阵；颜色越深表示正相关越强。 

par(op)

经验相关只用完整者，可能受选择影响；它用于提出结构候选，不应凭一张热图完成选择。应结合设计、时间间隔、参数数量、收敛和敏感性分析。

4 忽略相关性的代价

4.1 Naive OLS

naive_fit <- lm(
  sbp ~ program * time + I(time^2) + age_c + female,
  data = analysis_data
)
naive_term <- "programIntervention:time"
naive_table <- data.frame(
  方法 = "普通 OLS（错误地把记录视为独立）",
  估计 = coef(naive_fit)[naive_term],
  标准误 = coef(summary(naive_fit))[naive_term, "Std. Error"]
)
knitr::kable(naive_table, digits = 3, caption = "Naive OLS 的组别×时间项")
Naive OLS 的组别×时间项
方法 估计 标准误
programIntervention:time 普通 OLS（错误地把记录视为独立） -0.496 0.098

系数估计可能仍接近某个均值关系，但默认标准误把约 1642 行当作独立信息。重复测量共享参与者随机效应，因此这种精度声明不可信。

4.2 Independence-working 边际模型与 participant-cluster sandwich

下面从 OLS score 手写一个参与者聚类 sandwich。它适合展示“bread–meat–bread”逻辑，是 independence-working 的教学实现；它不是完整 GEE：没有迭代估计工作相关、small-sample CR2、缺失权重或稳健的有限样本校正。

cluster_vcov_cr1 <- function(model, cluster) {
  X <- model.matrix(model)
  u <- residuals(model)
  cluster <- droplevels(factor(cluster))
  groups <- levels(cluster)
  score <- t(vapply(groups, function(g) {
    idx <- which(cluster == g)
    as.vector(crossprod(X[idx, , drop = FALSE], u[idx]))
  }, numeric(ncol(X))))
  bread <- solve(crossprod(X))
  N <- nrow(X); p <- ncol(X); G <- length(groups)
  correction <- (G / (G - 1)) * ((N - 1) / (N - p))
  vc <- correction * bread %*% crossprod(score) %*% bread
  attr(vc, "clusters") <- G
  vc
}

vc_cr1 <- cluster_vcov_cr1(naive_fit, analysis_data$id)
cr1_se <- sqrt(diag(vc_cr1))[naive_term]
cr1_df <- attr(vc_cr1, "clusters") - 1
cr1_est <- coef(naive_fit)[naive_term]
cr1_result <- data.frame(
  方法 = c("OLS model-based", "Cluster sandwich CR1（教学版）"),
  估计 = cr1_est,
  标准误 = c(coef(summary(naive_fit))[naive_term, "Std. Error"], cr1_se),
  CI下限 = c(
    cr1_est - qt(0.975, df.residual(naive_fit)) * coef(summary(naive_fit))[naive_term, "Std. Error"],
    cr1_est - qt(0.975, cr1_df) * cr1_se
  ),
  CI上限 = c(
    cr1_est + qt(0.975, df.residual(naive_fit)) * coef(summary(naive_fit))[naive_term, "Std. Error"],
    cr1_est + qt(0.975, cr1_df) * cr1_se
  )
)
knitr::kable(cr1_result, digits = 3, caption = "独立 OLS 与 participant-cluster CR1 标准误")
独立 OLS 与 participant-cluster CR1 标准误
方法 估计 标准误 CI下限 CI上限
OLS model-based -0.496 0.098 -0.688 -0.303
Cluster sandwich CR1（教学版） -0.496 0.071 -0.635 -0.356

这里为教学展示使用 G−1G-1 个 cluster 自由度的 t 区间；这不等于 CR2/Satterthwaite，也不是所有设计下都最优。正式 GEE 通常使用专门软件，预先选择工作相关并报告稳健标准误；cluster 很少时必须采用适当的小样本方法。

Population-averaged 效应描述目标总体均值随组别和时间的变化；subject-specific 效应描述给定随机效应的个体条件关系。在线性 identity-link 模型中固定效应数值常相近；在 logistic 等非线性模型中二者一般不同，不能只换名称。

5 GLS：直接建模残差相关

5.1 CS 与 AR(1)

gls_formula <- sbp ~ program * time + I(time^2) + age_c + female
gls_cs <- nlme::gls(
  gls_formula, data = analysis_data, method = "REML", na.action = na.omit,
  correlation = nlme::corCompSymm(form = ~ 1 | id)
)
gls_ar1 <- nlme::gls(
  gls_formula, data = analysis_data, method = "REML", na.action = na.omit,
  correlation = nlme::corAR1(form = ~ visit_index | id)
)
gls_compare <- data.frame(
  结构 = c("CS", "AR(1)"),
  参数数 = c(attr(logLik(gls_cs), "df"), attr(logLik(gls_ar1), "df")),
  AIC = c(AIC(gls_cs), AIC(gls_ar1)),
  BIC = c(BIC(gls_cs), BIC(gls_ar1)),
  `组×月估计` = c(coef(gls_cs)[naive_term], coef(gls_ar1)[naive_term]),
  check.names = FALSE
)
knitr::kable(gls_compare, digits = 3, caption = "相同固定效应、REML 拟合下的 GLS 相关结构比较")
相同固定效应、REML 拟合下的 GLS 相关结构比较
结构 参数数 AIC BIC 组×月估计
CS 9 10095 10143 -0.459
AR(1) 9 9869 9917 -0.466

AIC/BIC 是相对拟合指标，不是相关结构为真的证明。比较固定效应相同的协方差结构可使用 ML 或 REML；比较不同固定效应时应使用 ML，确定固定效应后再用 REML 估计方差结构。AR(1) 的“一个 lag”在这里是一个计划访视（3 个月），不是任意实际天数。

5.2 MMRM：随访访视分类 + 基线调整 + 非结构化协方差

试验中一种常见 MMRM 把随访访视作为分类变量，加入组别×访视，调整同一人的基线 SBP，并让随访结局具有非结构化协方差。基线不是这一模型的响应；因此它与把五次访视都作为响应的 constrained longitudinal model 不同。

baseline_lookup <- long_data[
  long_data[["time"]] == 0, c("id", "sbp", "age_c", "female")
]
names(baseline_lookup)[2] <- "baseline_sbp"
baseline_lookup[["baseline_centered"]] <-
  baseline_lookup[["baseline_sbp"]] - mean(baseline_lookup[["baseline_sbp"]])
post_data <- merge(
  long_data[long_data[["time"]] > 0 & !is.na(long_data[["sbp"]]), ],
  baseline_lookup[c("id", "baseline_sbp", "baseline_centered")],
  by = "id", sort = FALSE
)
post_data <- post_data[order(as.integer(post_data[["id"]]), post_data[["time"]]), ]
post_data[["post_visit"]] <- factor(post_data[["time"]], levels = c(3, 6, 9, 12))
post_data[["post_index"]] <- match(post_data[["time"]], c(3, 6, 9, 12))

mmrm_fit <- nlme::gls(
  sbp ~ program * post_visit + baseline_centered + age_c + female,
  data = post_data, method = "REML", na.action = na.omit,
  correlation = nlme::corSymm(form = ~ post_index | id),
  weights = nlme::varIdent(form = ~ 1 | post_visit),
  control = nlme::glsControl(maxIter = 200, msMaxIter = 300)
)
mmrm_audit <- data.frame(
  项目 = c("有限 logLik", "有限系数", "apVar 可用", "相关参数数", "方差比参数数"),
  结果 = c(
    is.finite(as.numeric(logLik(mmrm_fit))), all(is.finite(coef(mmrm_fit))),
    is.matrix(mmrm_fit$apVar) && all(is.finite(mmrm_fit$apVar)),
    length(coef(mmrm_fit$modelStruct$corStruct, unconstrained = FALSE)),
    length(coef(mmrm_fit$modelStruct$varStruct, unconstrained = FALSE))
  )
)
knitr::kable(mmrm_audit, caption = "MMRM 数值与协方差参数审计")
MMRM 数值与协方差参数审计
项目 结果
有限 logLik 1
有限系数 1
apVar 可用 1
相关参数数 6
方差比参数数 3

这里将四个随访作为 factor，corSymm 估计 6 个相关，varIdent 允许每次随访残差方差不同。每人内的 post_index 是唯一整数。非结构化并不表示“没有假设”：它仍要求所有人共享同一个协方差矩阵，并依赖足够样本和稳定优化。正式监管分析还常规定有限样本自由度（如 Kenward–Roger）；nlme::gls() 的默认输出不是其替代。

5.3 标准化 MMRM 访视组差

mmrm_beta <- coef(mmrm_fit)
mmrm_V <- vcov(mmrm_fit)
mmrm_formula <- ~ program * post_visit + baseline_centered + age_c + female
mmrm_design <- function(program_level, visit_value) {
  base_people <- baseline_lookup
  nd <- data.frame(
    program = factor(program_level, levels = levels(long_data[["program"]])),
    post_visit = factor(visit_value, levels = c(3, 6, 9, 12)),
    baseline_centered = base_people[["baseline_centered"]],
    age_c = base_people[["age_c"]],
    female = base_people[["female"]]
  )
  X <- model.matrix(mmrm_formula, nd)
  colMeans(X[, names(mmrm_beta), drop = FALSE])
}
mmrm_contrasts <- do.call(rbind, lapply(c(3, 6, 9, 12), function(month) {
  L <- mmrm_design("Intervention", month) - mmrm_design("Usual care", month)
  est <- drop(L %*% mmrm_beta)
  se <- sqrt(drop(L %*% mmrm_V %*% L))
  data.frame(月 = month, `干预 − 常规` = est, 标准误 = se,
             `95% CI 下限` = est - 1.96 * se,
             `95% CI 上限` = est + 1.96 * se, check.names = FALSE)
}))
knitr::kable(mmrm_contrasts, digits = 2,
             caption = "标准化到共同基线协变量分布的随访组差（large-sample Wald CI）")
标准化到共同基线协变量分布的随访组差（large-sample Wald CI）
月 干预 − 常规 标准误 95% CI 下限 95% CI 上限
3 -1.41 0.44 -2.27 -0.54
6 -2.50 0.58 -3.64 -1.37
9 -4.39 0.63 -5.63 -3.15
12 -5.52 0.71 -6.90 -4.13

逐访视区间回答四个分别的问题；若四个访视都是 confirmatory，应预设多重性控制或联合检验。某次访视 P<0.05、另一次 P>0.05，并不能证明效应在两次访视之间突然“出现”或“消失”。

6 线性混合模型：个体轨迹异质性

6.1 随机截距、随机斜率与残差 AR(1)

lme_fit <- nlme::lme(
  fixed = sbp ~ program * time + I(time^2) + age_c + female,
  random = ~ time | id,
  correlation = nlme::corAR1(form = ~ visit_index | id),
  data = analysis_data, method = "REML", na.action = na.omit,
  control = nlme::lmeControl(opt = "optim", maxIter = 200, msMaxIter = 200)
)
lme_audit <- data.frame(
  项目 = c("有限 logLik", "固定效应均有限", "apVar 可用", "随机效应均有限"),
  通过 = c(
    is.finite(as.numeric(logLik(lme_fit))), all(is.finite(nlme::fixef(lme_fit))),
    is.matrix(lme_fit$apVar) && all(is.finite(lme_fit$apVar)),
    all(is.finite(as.matrix(nlme::ranef(lme_fit))))
  )
)
knitr::kable(lme_audit, caption = "混合模型拟合后的数值审计")
混合模型拟合后的数值审计
项目 通过
有限 logLik TRUE
固定效应均有限 TRUE
apVar 可用 TRUE
随机效应均有限 TRUE

仅仅没有报错不等于模型可靠。还要查看 apVar、方差是否贴近 0、相关是否贴近 ±1、固定效应与 BLUP 是否有限，以及替代起始值/结构是否给出相同结论。

REML 适合给定固定效应后的方差估计；若要比较“有无二次时间项”等不同固定效应，应把两个模型都用 method="ML" 重拟合，再用似然比、AIC 与科学合理性综合判断。不要用 REML likelihood 比较不同固定效应。

6.2 固定效应解释

lme_tab <- summary(lme_fit)$tTable
fixed_display <- data.frame(
  项 = rownames(lme_tab),
  估计 = lme_tab[, "Value"],
  标准误 = lme_tab[, "Std.Error"],
  自由度 = lme_tab[, "DF"],
  CI下限 = lme_tab[, "Value"] - qt(0.975, lme_tab[, "DF"]) * lme_tab[, "Std.Error"],
  CI上限 = lme_tab[, "Value"] + qt(0.975, lme_tab[, "DF"]) * lme_tab[, "Std.Error"],
  P值 = vapply(lme_tab[, "p-value"], format_p, character(1)),
  check.names = FALSE
)
knitr::kable(fixed_display, digits = 3, caption = "随机截距／斜率 + AR(1) 模型的固定效应")
随机截距／斜率 + AR(1) 模型的固定效应
项 估计 标准误 自由度 CI下限 CI上限 P值
(Intercept) (Intercept) 143.324 0.809 1279 141.738 144.911 <0.001
programIntervention programIntervention -0.090 0.928 356 -1.915 1.736 0.923
time time -0.257 0.082 1279 -0.417 -0.096 0.00173
I(time^2) I(time^2) 0.016 0.006 1279 0.004 0.027 0.00652
age_c age_c 1.185 0.408 356 0.383 1.988 0.00391
femaleYes femaleYes -2.463 0.837 356 -4.108 -0.817 0.00346
programIntervention:time programIntervention:time -0.465 0.064 1279 -0.590 -0.340 <0.001
  • 截距：55 岁、男性、常规照护者在第 0 月的条件平均 SBP；
  • programIntervention：第 0 月干预相对常规照护的调整组差，不是“总体干预效果”；
  • time：常规照护组在第 0 月的瞬时每月斜率；
  • programIntervention:time：干预组相对常规照护组的每月线性斜率差；
  • I(time^2)：两组共享的曲率；时点 tt 的斜率为相应线性部分加 2β̂t2t2\hat\beta_{t^2}t。

不要将 P>0.05 写成“无效”或“没有差异”。估计区间若同时容纳有益、无关紧要和有害值，正确结论是资料不够精确区分这些可能性。

6.3 用模型矩阵标准化绝对均值

系数表不直接给政策关心的各组各访视平均 SBP。下面把每个参与者的年龄和性别复制到每个 组 × 访视 情景，再平均模型矩阵；这样标准化到入组样本的共同协变量分布。

beta <- nlme::fixef(lme_fit)
V_beta <- vcov(lme_fit)
fixed_terms <- delete.response(terms(lme_fit))

standardized_row <- function(program_level, time_value) {
  nd <- data.frame(
    program = factor(program_level, levels = levels(long_data$program)),
    time = time_value,
    age_c = (age - 55) / 10,
    female = factor(female_num, levels = 0:1, labels = c("No", "Yes"))
  )
  X <- model.matrix(fixed_terms, nd)
  colMeans(X[, names(beta), drop = FALSE])
}

cell_grid <- expand.grid(
  program = levels(long_data$program), time = times,
  KEEP.OUT.ATTRS = FALSE, stringsAsFactors = FALSE
)
L_cells <- t(mapply(standardized_row, cell_grid$program, cell_grid$time,
                    SIMPLIFY = TRUE))
colnames(L_cells) <- names(beta)
cell_est <- as.vector(L_cells %*% beta)
cell_se <- sqrt(diag(L_cells %*% V_beta %*% t(L_cells)))
standardized_means <- data.frame(
  组别 = unname(zh_program[cell_grid$program]), 月 = cell_grid$time,
  调整均值 = cell_est,
  `95% CI 下限` = cell_est - 1.96 * cell_se,
  `95% CI 上限` = cell_est + 1.96 * cell_se,
  check.names = FALSE
)
knitr::kable(standardized_means, digits = 2,
             caption = "标准化到共同协变量分布的各组各访视 SBP 均值")
标准化到共同协变量分布的各组各访视 SBP 均值
组别 月 调整均值 95% CI 下限 95% CI 上限
常规照护 0 142 141 143
干预 0 142 140 143
常规照护 3 141 140 142
干预 3 140 139 141
常规照护 6 141 140 142
干预 6 138 137 139
常规照护 9 141 140 142
干预 9 137 135 138
常规照护 12 141 140 142
干预 12 135 134 137

这里使用固定效应协方差的 large-sample Wald 95% CI（±1.96 SE），不是个体预测区间。若目标总体的年龄/性别分布不同，应在目标总体上标准化，而非机械使用本样本。

6.4 12 月组差与 difference-in-changes

cell_key <- paste(cell_grid$program, cell_grid$time, sep = "@")
rownames(L_cells) <- cell_key
get_L <- function(g, t) L_cells[paste(g, t, sep = "@"), ]
L_u0 <- get_L("Usual care", 0); L_u12 <- get_L("Usual care", 12)
L_i0 <- get_L("Intervention", 0); L_i12 <- get_L("Intervention", 12)
contrast_list <- list(
  "12 月干预 − 常规" = L_i12 - L_u12,
  "常规组：12 月 − 基线" = L_u12 - L_u0,
  "干预组：12 月 − 基线" = L_i12 - L_i0,
  "变化差（干预 − 常规）" = (L_i12 - L_i0) - (L_u12 - L_u0)
)
contrast_table <- do.call(rbind, lapply(names(contrast_list), function(label) {
  L <- contrast_list[[label]]
  est <- sum(L * beta)
  se <- sqrt(drop(L %*% V_beta %*% L))
  data.frame(对比 = label, 估计 = est, 标准误 = se,
             `95% CI 下限` = est - 1.96 * se,
             `95% CI 上限` = est + 1.96 * se, check.names = FALSE)
}))
knitr::kable(contrast_table, digits = 2,
             caption = "预先定义的绝对 SBP 对比（large-sample Wald CI）")
预先定义的绝对 SBP 对比（large-sample Wald CI）
对比 估计 标准误 95% CI 下限 95% CI 上限
12 月干预 − 常规 -5.67 0.89 -7.42 -3.92
常规组：12 月 − 基线 -0.80 0.54 -1.87 0.26
干预组：12 月 − 基线 -6.38 0.54 -7.44 -5.33
变化差（干预 − 常规） -5.58 0.76 -7.08 -4.08

本模型允许随机组在基线有估计差异，所以 12 月组差与变化差不必完全相同。随机试验若预先采用 constrained longitudinal data analysis，可约束共同基线并提高效率，但必须在方案中明确。

7 预测与个体差异

7.1 Level 0 与 Level 1 预测

predict(..., level=0) 只用固定效应，回答同一协变量下的总体平均轨迹；level=1 加入已有参与者的经验 Bayes 随机效应（BLUP），回答给定其已观测资料后的条件轨迹。后者用于描述/预测已有个体，不能当作已知真值，也不能直接推广给新参与者。

completer_ids <- wide_sbp$id_num[complete.cases(wide_sbp)][1:6]
prediction_grid <- do.call(rbind, lapply(completer_ids, function(one_id) {
  d0 <- long_data[long_data$id_num == one_id, ][1, ]
  data.frame(id = factor(one_id, levels = levels(long_data$id)), id_num = one_id,
             time = times, visit_index = seq_along(times),
             program = factor(as.character(d0$program), levels = levels(long_data$program)),
             age_c = d0$age_c,
             female = factor(as.character(d0$female), levels = levels(long_data$female)))
}))
prediction_grid$level0 <- predict(lme_fit, newdata = prediction_grid, level = 0)
prediction_grid$level1 <- predict(lme_fit, newdata = prediction_grid, level = 1)
set_long_plot_font()
op <- par(mfrow = c(2, 3), mar = c(2.2, 2.2, 2.2, 0.8),
          oma = c(3.2, 3.5, 0.5, 0.5))
for (one_id in completer_ids) {
  obs <- analysis_data[analysis_data$id_num == one_id, ]
  pred <- prediction_grid[prediction_grid$id_num == one_id, ]
  plot(obs$time, obs$sbp, pch = 16, col = palette_long[["navy"]],
       xlab = "", ylab = "", main = paste("参与者", one_id),
       ylim = range(c(obs$sbp, pred$level0, pred$level1)))
  lines(pred$time, pred$level0, lty = 2, lwd = 2, col = palette_long[["gray"]])
  lines(pred$time, pred$level1, lwd = 2, col = palette_long[["orange"]])
}
mtext("月份", side = 1, outer = TRUE, line = 1.8)
mtext("SBP（mmHg）", side = 2, outer = TRUE, line = 2.0)
六个小面板展示个体收缩压观测点、总体预测虚线和个体条件预测实线。

六名完成者的观测值、总体固定效应轨迹（Level 0）和含 BLUP 的个体条件轨迹（Level 1）。 

par(op)

7.2 随机效应、时间变化的方差与相关

随机斜率使个体间方差随时间变化：

Var⁡(b0i+tb1i)=D00+2tD01+t2D11. \operatorname{Var}(b_{0i}+t b_{1i})=D_{00}+2tD_{01}+t^2D_{11}.

加上 AR(1) 残差后，两时点总协方差还包括 σ2ρ|j−k|\sigma^2\rho^{|j-k|}。因此“ICC”不再是一个常数；更清楚的做法是报告选定时点的随机轨迹方差占比及具体时点对的模型相关。

D_hat <- as.matrix(nlme::getVarCov(lme_fit, type = "random.effects"))
sigma2_hat <- sigma(lme_fit)^2
rho_hat <- as.numeric(coef(lme_fit$modelStruct$corStruct, unconstrained = FALSE))[1]
random_var <- function(t) drop(c(1, t) %*% D_hat %*% c(1, t))
total_var <- vapply(times, function(t) random_var(t) + sigma2_hat, numeric(1))
baseline_cov <- vapply(seq_along(times), function(j) {
  drop(c(1, 0) %*% D_hat %*% c(1, times[j])) + sigma2_hat * rho_hat^(j - 1)
}, numeric(1))
variance_table <- data.frame(
  月 = times,
  随机轨迹方差 = vapply(times, random_var, numeric(1)),
  残差方差 = sigma2_hat,
  随机轨迹方差占比 = vapply(times, random_var, numeric(1)) / total_var,
  与基线模型相关 = baseline_cov / sqrt(total_var[1] * total_var)
)
knitr::kable(variance_table, digits = 3,
             caption = "随机斜率模型下随时间变化的方差分解与基线相关")
随机斜率模型下随时间变化的方差分解与基线相关
月 随机轨迹方差 残差方差 随机轨迹方差占比 与基线模型相关
0 60.7 18.4 0.768 1.000
3 55.6 18.4 0.752 0.882
6 52.2 18.4 0.740 0.802
9 50.5 18.4 0.733 0.738
12 50.5 18.4 0.733 0.680

8 模型诊断与稳健性

8.1 残差和随机效应检查

std_resid <- residuals(lme_fit, type = "normalized")
fit_cond <- fitted(lme_fit, level = 1)
re_hat <- nlme::ranef(lme_fit)
set_long_plot_font()
op <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(fit_cond, std_resid, pch = 16, cex = 0.5,
     col = grDevices::adjustcolor(palette_long[["blue"]], 0.4),
     xlab = "条件拟合值", ylab = "标准化残差", main = "残差 vs 拟合值")
abline(h = 0, lty = 2, col = palette_long[["gray"]])
plot(analysis_data$time, std_resid, pch = 16, cex = 0.5,
     col = grDevices::adjustcolor(palette_long[["teal"]], 0.35), xlab = "月份", ylab = "标准化残差",
     main = "残差 vs 时间")
abline(h = 0, lty = 2, col = palette_long[["gray"]])
qqnorm(std_resid, pch = 16, cex = 0.45,
       col = grDevices::adjustcolor(palette_long[["navy"]], 0.4),
       main = "残差 Q-Q 图", xlab = "理论正态分位数", ylab = "样本分位数")
qqline(std_resid, col = palette_long[["vermillion"]], lwd = 2)
plot(re_hat[, 1], re_hat[, 2], pch = 16, cex = 0.55,
     col = grDevices::adjustcolor(palette_long[["orange"]], 0.45),
     xlab = "随机截距 BLUP", ylab = "随机斜率 BLUP", main = "随机效应")
abline(h = 0, v = 0, lty = 3, col = palette_long[["gray"]])
四个诊断面板依次显示残差对拟合值、残差对月份、残差正态分位图和随机截距斜率散点。

混合模型的标准化条件残差与随机效应诊断；诊断图用于发现模型偏离而非机械通过测试。 

par(op)

还应检查：异常个体删除敏感性、按组/访视的残差方差、实际时间间隔、随机结构简化、CS/AR(1)/非结构化结果一致性、非线性、结局变换，以及是否有少数参与者主导随机斜率。正态性主要针对条件残差和随机效应分布，不要求原始 SBP 完全正态。

9 缺失数据：模型能做什么、不能做什么

9.1 MCAR、MAR 与 MNAR

机制 定义（相对于分析中可用资料） 主要含义
MCAR 缺失与已观测和未观测结局都无关 完整病例可无偏但低效；仍应报告
MAR 给定已观测历史/协变量后，与未观测结局无关 正确似然模型或适当 MI 可有效
MNAR 即使给定已观测资料，缺失仍依赖未观测结局 需不可检验假设和敏感性分析

本页退出概率只依赖上一时点已观测 SBP 与随机组；仍 active 的人上一值必定已观测，退出后不再抽签。因此 DGM 是给定该历史的 MAR。lme()/gls() 的 observed-data likelihood 在模型和 MAR 正确时可使用所有可用结局；它们不会“自动处理 MNAR”。

9.2 完整病例陷阱

cc_ids <- wide_sbp$id_num[complete.cases(wide_sbp)]
cc_data <- analysis_data[analysis_data$id_num %in% cc_ids, ]
cc_fit <- lm(sbp ~ program * time + I(time^2) + age_c + female, data = cc_data)
missing_compare <- data.frame(
  分析 = c("所有观测 + LMM likelihood", "五次均完整 + OLS"),
  参与者数 = c(length(unique(analysis_data$id)), length(unique(cc_data$id))),
  `组×月估计` = c(nlme::fixef(lme_fit)[naive_term], coef(cc_fit)[naive_term]),
  标准误 = c(summary(lme_fit)$tTable[naive_term, "Std.Error"],
             coef(summary(cc_fit))[naive_term, "Std. Error"])
)
knitr::kable(missing_compare, digits = 3,
             caption = "可用资料似然分析与完整五访视分析的比较")
可用资料似然分析与完整五访视分析的比较
分析 参与者数 组.月估计 标准误
所有观测 + LMM likelihood 360 -0.465 0.064
五次均完整 + OLS 300 -0.437 0.101

完整病例把 estimand 改成“会完成全部访视者”的轨迹；如果完成与既往 SBP 有关，这群人不是随机子样本。主要分析外还应做：

  • MI：插补模型包含结局历史、随机组、时间、强预测因子及分析模型的交互/非线性；反映多层结构；
  • IPW：建模每次继续被观测概率，使用稳定权重；检查极端权重和 positivity；
  • Pattern-mixture / delta adjustment / tipping point：系统改变未观测结局相对 MAR 插补值；
  • Selection 或 joint model：用明确的退出/结局共享机制；结果依赖额外假设；
  • 死亡等 intercurrent event：不能一概作为普通缺失，应先定义 treatment-policy、composite、while-alive 等 estimand。

MI、IPW 和 joint model 不是可互换的按钮。方法、变量、时间顺序、权重截尾和 MNAR sensitivity 参数应在看结果前规定。

10 Time-varying covariate

10.1 Within–between 分解

假设每次已完成访视还记录活动分钟。下面的教学变量只在结局实际观测行构造，不使用退出后的 full_sbp；真实研究应直接使用按时间测得的活动资料。

tv_data <- analysis_data
tv_data$activity <- pmax(
  0, 52 + 0.9 * tv_data$time + 4 * tv_data$program_num - 0.15 * (tv_data$sbp - 140)
)
tv_data$activity_between <- ave(tv_data$activity, tv_data$id, FUN = mean)
tv_data$activity_within <- tv_data$activity - tv_data$activity_between
tv_fit <- nlme::lme(
  sbp ~ program * time + I(time^2) + age_c + female +
    activity_within + activity_between,
  random = ~ time | id, data = tv_data, method = "REML",
  control = nlme::lmeControl(opt = "optim")
)
tv_terms <- summary(tv_fit)$tTable[c("activity_within", "activity_between"), , drop = FALSE]
knitr::kable(data.frame(效应 = c("个体内偏离", "个体间平均差"),
                        估计 = tv_terms[, "Value"], 标准误 = tv_terms[, "Std.Error"]),
             digits = 3, caption = "Time-varying activity 的 within–between 分解示例")
Time-varying activity 的 within–between 分解示例
效应 估计 标准误
activity_within 个体内偏离 -6.67 0
activity_between 个体间平均差 -6.67 0

activity_within 比较同一个人在“比自己通常水平更活跃”的访视，activity_between 比较长期平均活动不同的人。它们不是同一个问题。这里活动由同访视 SBP 构造，故存在同时性/反向因果，只用于语法教学。

尤其在随机试验中，活动可能受干预影响。把 post-randomization activity 加入主要模型会偏离随机分配总效应这个 estimand；普通回归调整通常不能识别直接效应，并可能产生 collider bias。若目标是中介效应，应使用明确的纵向因果框架，而非把它当“更充分调整”。

11 非连续纵向结局

11.1 二元结局：GEE 与 logistic GLMM

对反复“血压是否控制”，population-averaged logistic GEE 与带个体随机截距的 logistic GLMM 估计不同尺度的 odds ratio。下面代码只示范接口，不在本页生成额外随机结局；运行前检查 lme4。

if (!requireNamespace("lme4", quietly = TRUE)) {
  stop("此扩展示例需要 lme4 包。")
}
binary_data <- analysis_data
binary_data$controlled <- as.integer(binary_data$sbp < 130)
binary_glmm <- lme4::glmer(
  controlled ~ program * time + I(time^2) + age_c + female + (1 | id),
  family = binomial(), data = binary_data
)
# 固定效应 exp(beta) 是 subject-specific 条件 odds ratio；
# 应另用标准化预测报告各组绝对概率和概率差。

完整 GEE 可用专门软件指定 id、工作相关与稳健标准误。失访与结局有关时，普通 GEE 的完全随机缺失条件可能不足，需要加权 GEE 或多重插补，并明确假设。

11.2 计数结局：Poisson/负二项 GLMM

if (!requireNamespace("lme4", quietly = TRUE)) {
  stop("此扩展示例需要 lme4 包。")
}
set.seed(20261007)
count_data <- data.frame(
  id = factor(rep(seq_len(200), each = 4)),
  time = rep(c(0, 1, 2, 3), times = 200),
  program = factor(rep(rep(c("Usual care", "Intervention"), each = 100), each = 4)),
  person_months = runif(800, 0.7, 1.0)
)
# 数据契约：每行是一段 person-period；visits 是该段计数，person_months > 0 是暴露时长。
count_random_intercept <- rep(rnorm(200, 0, 0.45), each = 4)
count_mean <- with(
  count_data,
  person_months * exp(-0.2 - 0.18 * (program == "Intervention") * time +
                        0.08 * time + count_random_intercept)
)
count_data$visits <- rpois(nrow(count_data), count_mean)
count_glmm <- lme4::glmer(
  visits ~ program * time + offset(log(person_months)) + (1 | id),
  family = poisson(), data = count_data
)
# 检查过度离散、零膨胀、暴露时间与 dropout；必要时考虑负二项或其他模型。

计数模型的指数化系数是条件 rate ratio；offset 系数固定为 1。无论二元或计数，都应从模型给出绝对风险/率、明确边际或条件尺度，并让相关层级与设计一致。

12 样本量与功效

纵向研究的功效取决于主要对比、访视数/间隔、个体间异质性、相关结构、残差方差、失访和分析模型。不能先看本研究 P 值再报告“observed power”；它主要是 P 值的重新表达。

planned_per_arm <- seq(80, 240, by = 20)
dropout_scenarios <- c(0.10, 0.20, 0.30)
set_long_plot_font()
op <- par(mar = c(4.4, 4.4, 2.5, 1), xpd = NA)
matplot(planned_per_arm,
        outer(planned_per_arm, 1 - dropout_scenarios), type = "l", lwd = 2, lty = 1,
        col = c(palette_long[["green"]], palette_long[["orange"]], palette_long[["vermillion"]]),
        xlab = "每组计划入组人数", ylab = "预期 12 月保留人数",
        ylim = c(0, max(planned_per_arm)))
legend("top", inset = -0.12, horiz = TRUE, bty = "n",
       legend = paste0("失访 ", 100 * dropout_scenarios, "%"),
       col = c(palette_long[["green"]], palette_long[["orange"]], palette_long[["vermillion"]]), lwd = 2)
三条曲线显示每组计划入组数增加时，在10%、20%和30%失访率下预期保留人数。

在不同假设失访率下，每组计划入组数与预期 12 月保留人数的关系；这是保留率情景图，不是完整功效计算。 

par(op)

正式规划应围绕 12 月组差、斜率差或曲线下面积等一个主要 estimand，使用外部/先导资料设定协方差，并模拟完整流程：生成随机效应和残差、应用访视窗口和失访、拟合计划模型、记录 CI 覆盖和拒绝比例。还要改变效应、相关、MNAR 偏移和 cluster 数做情景分析；模拟次数要报告 Monte Carlo 误差。

13 从方案到报告

13.1 可审核分析工作流

  1. 写清目标总体、组别、结局、访视、效应尺度和 intercurrent-event 策略。
  2. 确认随机化/抽样单位、重复单位和 long 数据唯一键。
  3. 冻结时间函数、基线处理、主要对比与目标标准化总体。
  4. 画个体轨迹、各组绝对均值、保留率和访视模式，不按 P 值筛变量。
  5. 依据设计提出 CS、AR(1)、非结构化或随机斜率候选；限制复杂度。
  6. 用 ML 比较不同固定效应；固定结构后用 REML 估计连续结局 LMM。
  7. 报告绝对标准化均值、预设对比、CI 和 clinically meaningful threshold。
  8. 检查残差、随机效应、收敛、边界方差、异常个体和替代结构。
  9. 陈述缺失假设；按预案实施 MI/IPW/MNAR sensitivity，而非只做完整病例。
  10. 保存代码、种子、软件版本、模型公式、协方差参数和所有分析偏离。

13.2 报告模板

共随机分配 360 名参与者，在基线及第 3、6、9、12 月计划测量 SBP；12 月保留率为 …。主要 estimand 为入组总体中干预相对常规照护的 12 月平均 SBP 差。使用包含随机组、连续时间、二次时间、随机组×时间、预设基线协变量、个体随机截距/斜率及访视级 AR(1) 残差的线性混合模型，并在 MAR（给定已观测历史）下用可用资料 REML 估计。标准化 12 月均值为 … 和 … mmHg，组差为 … mmHg（95% CI …）。替代 CS、MMRM 与完整病例结果为 …。诊断提示 …。结论依赖模型、MAR、测量和实施假设；MNAR tipping-point 分析显示 …。

13.3 常见错误与修正

常见错误 问题 修正
把所有记录视为独立 标准误忽略人内相关 GEE/GLS/LMM 或合适 cluster-robust 方法
只比较每组各自 P 值 “一组显著、另一组不显著”不等于组差 直接估计 program × time 或预设对比
把 program 主项叫全程效果 有交互时它是 time=0 差 报告各时点标准化均值与明确对比
用 REML 比较不同固定效应 REML likelihood 基于不同变换 固定效应比较用 ML；选定后 REML
先试很多协方差再挑最小 P 结果驱动选择低估不确定性 预设少量候选，报告敏感性
删除退出者 改变目标人群且可能偏倚 可用资料似然 + 缺失敏感性分析
说 LMM 自动处理 MNAR 似然通常依赖 MAR 明确 MAR，增加 MNAR sensitivity
调整干预后的协变量 可能改变 estimand/引入偏差 总效应模型不机械调整；中介另定义
P>0.05 宣称无效 可能只是区间宽 报估计、CI 和重要效应阈值
把 BLUP 当真实个体效应 BLUP 被收缩且有估计误差 标为条件预测，给预测不确定性

14 知识检查与练习

快速自测

  1. 五个年份分别抽取不同居民，为什么不能估计居民个体内变化？
  2. 在 program * time 中，programIntervention 代表哪个时点的组差？
  3. AR(1) 的 ρ=0.6\rho=0.6 时，相隔两个计划访视的残差相关是多少？
  4. 为什么“LMM 使用所有可用资料”不等于“对 MNAR 无偏”？
  5. Level 0 和 Level 1 预测分别适合什么对象？
  6. 活动量若受干预影响，为何不能随意加入总效应模型？
查看答案
  1. 同一个人没有重复测量；只能估计各年份总体构成和均值的变化，个体变化不可识别。
  2. time=0 的干预相对常规照护差；交互项才是线性斜率差。
  3. 0.62=0.360.6^2=0.36，前提是“两个 lag”按相同计划间隔定义。
  4. Observed-data likelihood 的无偏性仍要求模型正确和给定已观测资料的 MAR；MNAR 涉及未观测值，需要额外不可检验假设。
  5. Level 0 是固定效应总体轨迹，适合目标总体或新个体的均值；Level 1 加已有个体 BLUP，适合已有观测者的条件轨迹。
  6. 它是 post-randomization 变量，可能是中介或 collider；调整会改变 estimand，不能再称随机分配的总效应。

实践练习

  1. 将 time 中心化到 12 月并重拟合，验证 program 主项等于 12 月组差，而拟合值不变。
  2. 用 MMRM 的模型矩阵估计五个访视的组差；与二次 LMM 对比，讨论形状假设与精度。
  3. 把 AR(1) 改为 CS，比较 AIC、残差图和 12 月对比；不要只比较 P 值。
  4. 只保留完成五访视者，说明分析人群、估计和标准误如何改变，以及为何不能据此判断 MAR/MNAR。
  5. 为“退出后未观测 SBP 比 MAR 预测高 2、4、6 mmHg”设计 delta-adjustment sensitivity table。
  6. 设计一项以 12 月组差为主要 estimand 的模拟功效研究，至少改变失访、随机斜率 SD 和 AR(1) 相关。
练习讨论要点
  1. 使用 time12=time-12；交互斜率和预测不变，截距/组主项改变解释。
  2. MMRM 为每个访视单独估计，不强迫二次曲线；参数更多。比较绝对估计和 CI，而非胜负标签。
  3. 固定效应相同可比较 ML/REML 信息准则；同时检查科学合理性、稳定性和主要对比敏感度。
  4. 完整者由退出过程选择，结果针对完成者；观测资料无法检验未观测结局依赖关系。
  5. 对每个 delta 重做插补/估计并画组差；找结论跨越临床界值或零点的 tipping point。
  6. 每次生成完整轨迹、施加失访、拟合预设模型；功效为预设检验拒绝比例，并报告 Monte Carlo SE、偏倚和 CI 覆盖。

15 速查表与提交前清单

速查表

目标 R 语法 解释提醒
Long/wide 转换 reshape(..., direction=) 转换前验证 id × time 唯一
描述保留率 aggregate(observed ~ time, mean) 按组和退出原因继续分层
Independence working lm() + cluster sandwich 教学 CR1 不等于完整 GEE/CR2
GLS–CS gls(..., corCompSymm()) 所有时点对同一相关
GLS–AR(1) gls(..., corAR1()) lag 定义须对应时间设计
MMRM gls(..., corSymm(), varIdent()) 参数多；有限样本 df 另处理
随机轨迹 lme(..., random=~time|id) 固定效应条件于随机效应
总体预测 predict(fit, level=0) 不含个体 BLUP
个体条件预测 predict(fit, level=1) 仅对已有个体、存在收缩
固定效应比较 method="ML" 不同固定效应勿比 REML likelihood
可用资料似然 na.action=na.omit 正确模型 + MAR，不涵盖 MNAR
Within–between x-xbar 与 xbar 两个不同 estimand

提交前清单

## 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   xfun_0.60       lattice_0.22-9 
##  [6] cachem_1.1.0    knitr_1.51      htmltools_0.5.9 rmarkdown_2.31  stats4_4.6.1   
## [11] lifecycle_1.0.5 cli_3.6.6       grid_4.6.1      sass_0.4.10     jquerylib_0.1.4
## [16] compiler_4.6.1  tools_4.6.1     nlme_3.1-169    evaluate_1.0.5  bslib_0.12.0   
## [21] yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0