The V Lab
适用对象医学、公共卫生、社会学、心理学与健康数据科学学习者
学习时长约 180–240 分钟
先修要求线性／逻辑回归、置信区间、方差与基本 R 语法

关于本教程的数据 患者、诊所、医院、学生、学校、心理治疗师和随访记录均由固定随机种子模拟，不含真实个人信息。模拟机制用于验证代码和解释，不构成任何真实治疗、机构或社会政策证据。

如何使用本教程

多层模型不应从一条复杂的 lmer() 公式开始。推荐按以下顺序学习：

确认分析单位与依赖结构 → 写出目标关系 → 分解方差 → 区分层内与层间效应 → 选择随机结构 → 拟合与诊断 → 预测、敏感性分析和报告

本页用四套可复现模拟数据贯穿讲解：患者嵌套于诊所的连续结局、学生嵌套于学校的随机斜率模型、治疗时点嵌套于个体再嵌套于治疗师的三层模型，以及患者嵌套于医院的二分类结局。所有主要代码只依赖 R 自带功能、lme4 和页面渲染所需的 knitr。

学习目标

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

  • 根据研究设计识别观测、个体和群组层级，以及嵌套、交叉分类和重复测量；
  • 解释固定效应、随机效应、方差成分与部分汇聚，而不把“随机”理解为“不重要”；
  • 由随机截距模型计算并正确解释 ICC；
  • 区分完全汇聚、不汇聚与部分汇聚；
  • 使用群组均值分解区分 within 与 between 关系；
  • 解释随机斜率、截距–斜率协方差和跨层交互；
  • 为连续、二分类和纵向结局选择相应的 LMM 或 GLMM；
  • 区分 cluster-specific 条件效应与总体边际效应；
  • 正确选择 ML 与 REML，并识别收敛、奇异拟合和影响群组；
  • 说明 mixed model 对缺失、序列相关和因果识别的边界；
  • 使用 lme4::lmer()、lme4::glmer()、lme4::VarCorr() 与预测接口；
  • 写出包含样本层级、固定效应、随机结构、诊断与局限的可审核报告。

1 为什么普通回归可能不够

1.1 独立性由数据生成过程决定

同一诊所的患者可能共享转诊标准、医生、设备和地区环境；同一学校的学生共享课程和资源；同一个人的日记时点共享稳定人格和经历。即使每一行看起来都是一条完整记录，这些共享因素也会使误差相关。

若忽略这种依赖：

  • 标准误和置信区间可能错误；
  • 不同层级的关系会被混在一起；
  • 群组间异质性被压进一个残差项；
  • 预测新个体与预测已有个体会被混为一谈；
  • 样本量看似很大，但高层变量的独立信息仍可能很少。
design_map <- data.frame(
  领域 = c("医学：单次结局", "医学：重复随访", "社会学", "心理学", "医学：二分类"),
  `Level 1` = c("患者", "随访时点", "学生", "治疗时点", "患者"),
  `Level 2` = c("诊所", "患者", "学校", "个体", "医院"),
  `Level 3／交叉结构` = c("—", "医院", "社区可与学校交叉", "治疗师", "—"),
  主要问题 = c(
    "诊所内患者是否相关？",
    "恢复轨迹是否因患者与医院不同？",
    "个体 SES 与学校构成效应是否不同？",
    "当天压力与长期压力是否不同？",
    "医院差异与干预的条件 OR 是多少？"
  ),
  check.names = FALSE
)

knitr::kable(design_map, caption = "三个应用领域中的层级、单位与问题")
三个应用领域中的层级、单位与问题
领域 Level 1 Level 2 Level 3／交叉结构 主要问题
医学：单次结局 患者 诊所 — 诊所内患者是否相关？
医学：重复随访 随访时点 患者 医院 恢复轨迹是否因患者与医院不同？
社会学 学生 学校 社区可与学校交叉 个体 SES 与学校构成效应是否不同？
心理学 治疗时点 个体 治疗师 当天压力与长期压力是否不同？
医学：二分类 患者 医院 — 医院差异与干预的条件 OR 是多少？

有 ID 不等于自动需要随机截距 “层级”是相对于研究问题和依赖机制定义的，不是变量固有属性。必须说明哪些记录共享什么来源，以及结果要推广到个体、群组还是两者；不能只因为数据含有 hospital_id 就机械添加 (1 | hospital_id)。

1.2 嵌套、交叉与多重成员关系

  • 嵌套：每名患者只属于一家诊所；每个时点只属于一个人。
  • 交叉分类：学生同时属于学校与居住社区，而学校并不严格属于单一社区。
  • 多重成员关系：一名患者在疗程中由多名治疗师共同服务。
  • 不平衡群组：不同医院患者数不同；这在真实数据中很常见。
  • 信息性群组大小：群组规模本身与结局风险相关，普通模型可能需要额外处理。
检验你的理解：5000 名患者来自 8 家医院，医院层样本量是多少？ 对医院层暴露或政策效应而言，独立信息主要来自 8 家医院，而不是 5000 名患者。大量患者能改善医院内均值估计，却不能创造更多独立医院。

2 两层随机截距模型

2.1 固定关系与群组偏离

令个体 (i) 嵌套于群组 (j)：

Yij=β0+β1Xij+β2Zj+u0j+εij, Y_{ij}=\beta_0+\beta_1X_{ij}+\beta_2Z_j+u_{0j}+\varepsilon_{ij}, u0j∼N(0,τ00),εij∼N(0,σ2). u_{0j}\sim N(0,\tau_{00}),\qquad \varepsilon_{ij}\sim N(0,\sigma^2).

XijX_{ij} 是个体层变量，ZjZ_j 是群组层变量。固定效应 β\beta 描述所有群组共享的平均条件关系；随机截距 u0ju_{0j} 描述第 jj 个群组相对总体截距的偏离，其分布方差为 τ00\tau_{00}。

“固定”和“随机”不是日常语言 固定效应不表示变量不会变化；随机效应也不表示变量不重要或群组一定由简单随机抽样获得。这里的区别是：模型直接估计共同系数 β\beta，同时用一个分布描述许多群组偏离 uju_j。

2.2 ICC 是相关结构，不是因果比例

空随机截距模型中的组内相关系数为：

ICC=ρ=τ00τ00+σ2. ICC=\rho=\frac{\tau_{00}}{\tau_{00}+\sigma^2}.

它等于同一群组随机抽取两名个体的模型相关性，也可描述空模型总方差中群组间成分的比例。它不能解释为“群组造成了该百分比的结局”。调整协变量后计算的 conditional ICC 回答的是另一个方差分解问题。

低 ICC 也不自动允许忽略聚类：影响还取决于平均群组大小、群组数、暴露所在层级和不平衡程度。

2.3 完全汇聚、不汇聚与部分汇聚

策略 做法 优点 局限
完全汇聚 忽略群组差异 简单 把相关性和异质性压入残差
不汇聚 每个群组单独估计 保留群组差异 小群组极不稳定，参数很多
部分汇聚 群组估计向总体适度收缩 稳定且保留异质性 依赖 random-effect 分布与模型

简单空模型下，群组偏离可近似表示为：

ûj≈λj(Y‾j−β̂0),λj=τ00τ00+σ2/nj. \widehat u_j\approx \lambda_j(\bar Y_j-\widehat\beta_0),\qquad \lambda_j=\frac{\tau_{00}}{\tau_{00}+\sigma^2/n_j}.

小群组信息较少，λj\lambda_j 较小，收缩更强；大群组更多保留自身信息。部分汇聚不是“把机构差异平均掉”，而是在群组数据与总体分布间进行有原则的权衡。

3 医学应用一：患者嵌套于诊所

3.1 模拟队列与时间零点

设 48 家诊所各纳入 22–36 名患者。干预在诊所层分配，结局为 6 个月收缩压；患者治疗前记录年龄和基线收缩压。由于治疗在诊所层变化，治疗效应的高层独立信息来自诊所数。

set.seed(20260821)
J <- 48L
n_j <- sample(22:36, J, replace = TRUE)
clinic <- sprintf("C%02d", seq_len(J))
treat_j <- sample(rep(0:1, each = J / 2))
u0 <- rnorm(J, 0, 4.5)

med <- do.call(rbind, lapply(seq_len(J), function(j) {
  n <- n_j[j]
  age <- pmin(pmax(rnorm(n, 60, 11), 30), 85)
  baseline <- rnorm(n, 138 + 0.08 * (age - 60), 12)
  sbp_6m <- 132 - 4.2 * treat_j[j] +
    0.55 * (baseline - 138) + 0.12 * (age - 60) +
    u0[j] + rnorm(n, 0, 9)
  data.frame(
    clinic = clinic[j],
    treatment = treat_j[j],
    age_c10 = (age - 60) / 10,
    baseline_c10 = (baseline - 138) / 10,
    sbp_6m = sbp_6m
  )
}))
med$treatment <- factor(
  med$treatment,
  0:1,
  c("Usual care", "Intervention")
)

med_audit <- data.frame(
  患者数 = nrow(med),
  诊所数 = length(unique(med$clinic)),
  最小诊所人数 = min(table(med$clinic)),
  中位诊所人数 = median(table(med$clinic)),
  最大诊所人数 = max(table(med$clinic)),
  check.names = FALSE
)

knitr::kable(med_audit, caption = "患者—诊所模拟队列审计")
患者—诊所模拟队列审计
患者数 诊所数 最小诊所人数 中位诊所人数 最大诊所人数
1425 48 22 29.5 36

3.2 空模型、ICC 与部分汇聚

m_med_empty <- lme4::lmer(
  sbp_6m ~ 1 + (1 | clinic),
  data = med,
  REML = TRUE
)

med_empty_icc <- icc_random_intercept(m_med_empty, "clinic")
med_empty_vc <- as.data.frame(lme4::VarCorr(m_med_empty))

med_empty_variance <- data.frame(
  成分 = c("诊所随机截距", "患者层残差"),
  方差 = c(
    med_empty_vc$vcov[
      med_empty_vc$grp == "clinic" & is.na(med_empty_vc$var2)
    ],
    stats::sigma(m_med_empty)^2
  ),
  标准差 = sqrt(c(
    med_empty_vc$vcov[
      med_empty_vc$grp == "clinic" & is.na(med_empty_vc$var2)
    ],
    stats::sigma(m_med_empty)^2
  )),
  check.names = FALSE
)

knitr::kable(
  med_empty_variance,
  digits = 3,
  caption = "医学空随机截距模型的方差成分"
)
医学空随机截距模型的方差成分
成分 方差 标准差
诊所随机截距 27.1 5.21
患者层残差 128.6 11.34

空模型 ICC 为 0.174：同一诊所两名随机患者的 6 个月收缩压在模型尺度上具有约 17.4% 的相关性。它不表示诊所因果造成了 17.4% 的血压。

clinic_raw <- aggregate(sbp_6m ~ clinic, data = med, FUN = mean)
clinic_raw$n <- as.numeric(table(med$clinic)[clinic_raw$clinic])
clinic_re <- lme4::ranef(m_med_empty)$clinic
clinic_partial <- data.frame(
  clinic = rownames(clinic_re),
  partial = unname(lme4::fixef(m_med_empty)[1] + clinic_re[, 1])
)
clinic_compare <- merge(clinic_raw, clinic_partial, by = "clinic")
clinic_compare <- clinic_compare[order(clinic_compare$sbp_6m), ]
plot_y <- seq_len(nrow(clinic_compare))

plot(
  clinic_compare$sbp_6m,
  plot_y,
  pch = 16,
  col = palette_ml["orange"],
  xlim = range(c(clinic_compare$sbp_6m, clinic_compare$partial)),
  yaxt = "n",
  xlab = "6 个月平均收缩压",
  ylab = "按原始均值排序的诊所",
  main = "原始均值与部分汇聚",
  las = 1
)
segments(
  clinic_compare$sbp_6m,
  plot_y,
  clinic_compare$partial,
  plot_y,
  col = palette_ml["gray"]
)
points(
  clinic_compare$partial,
  plot_y,
  pch = 17,
  col = palette_ml["teal"]
)
abline(v = lme4::fixef(m_med_empty)[1], lty = 2, col = palette_ml["navy"])
legend(
  "bottomright",
  legend = c("原始诊所均值", "部分汇聚均值", "总体截距"),
  col = c(palette_ml["orange"], palette_ml["teal"], palette_ml["navy"]),
  pch = c(16, 17, NA),
  lty = c(NA, NA, 2),
  bty = "n"
)
四十八家诊所按原始平均收缩压排序。橙色点表示原始均值，蓝绿色点表示部分汇聚估计，连接线大多指向总体均值，小诊所的移动通常更明显。

诊所原始均值与空随机截距模型的部分汇聚均值。线段连接同一家诊所的两个估计；小样本诊所通常向总体均值收缩更多。

条件模态常被称为 BLUP，但不是没有误差的“诊所真值”。将其用于机构排名、处罚或资源配置时，必须同时呈现不确定性、病例结构和收缩程度。

3.3 调整后诊所随机截距模型

m_med <- lme4::lmer(
  sbp_6m ~ treatment + baseline_c10 + age_c10 + (1 | clinic),
  data = med,
  REML = TRUE
)

med_labels <- c(
  `(Intercept)` = "常规管理、基线协变量为中心值",
  treatmentIntervention = "干预 vs 常规管理",
  baseline_c10 = "基线收缩压每增加 10 mmHg",
  age_c10 = "年龄每增加 10 岁"
)
med_fixed <- lmm_fixed_table(m_med, med_labels)
med_adjusted_icc <- icc_random_intercept(m_med, "clinic")

knitr::kable(
  med_fixed,
  digits = 3,
  caption = "患者—诊所调整后线性混合模型固定效应"
)
患者—诊所调整后线性混合模型固定效应
变量 估计值 标准误 95% CI 下限 95% CI 上限 t 值
(Intercept) 常规管理、基线协变量为中心值 132.219 1.031 130.199 134.24 128.27
treatmentIntervention 干预 vs 常规管理 -4.113 1.457 -6.969 -1.26 -2.82
baseline_c10 基线收缩压每增加 10 mmHg 5.568 0.217 5.143 5.99 25.71
age_c10 年龄每增加 10 岁 0.573 0.228 0.127 1.02 2.52

调整后干预组平均 6 个月收缩压比常规管理低 4.11 mmHg（Wald 95% CI：-6.97 至 -1.26）。这是给定基线收缩压、年龄和诊所随机截距后的条件均值差。由于干预在诊所层分配，因果解释还需要诊所层可交换性、足够重叠和清晰的分配机制。

调整后 ICC 为 0.207，高于空模型的 0.174 并不矛盾：移除部分患者层可解释方差后，剩余总方差的构成可以改变。ICC 不是必须随协变量增加而下降的“模型质量分数”。

med_residuals <- residuals(m_med, type = "pearson")
med_fitted <- fitted(m_med)
old_par <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 3, 1))

plot(
  med_fitted,
  med_residuals,
  pch = 16,
  cex = 0.55,
  col = grDevices::adjustcolor(palette_ml["teal"], alpha.f = 0.45),
  xlab = "条件拟合值",
  ylab = "Pearson 残差",
  main = "残差与拟合值",
  las = 1
)
abline(h = 0, lty = 2, col = palette_ml["vermillion"])

stats::qqnorm(
  as.numeric(scale(med_residuals)),
  pch = 16,
  cex = 0.55,
  col = palette_ml["blue"],
  main = "残差 Q–Q 图"
)
stats::qqline(as.numeric(scale(med_residuals)), col = palette_ml["vermillion"])
两个诊断面板。左侧散点围绕零水平线分布，右侧标准化残差点大致沿正态参考直线排列，并在尾部存在随机偏离。

患者—诊所线性混合模型的条件残差诊断。左图比较拟合值与残差，右图为标准化残差正态 Q–Q 图；诊断用于寻找非线性、异方差和尾部偏离。 

par(old_par)

残差图不是唯一诊断。还应检查诊所层影响点、诊所大小、random-effect 分布、基线收缩压的函数形式，以及是否存在诊所特异残差方差。

4 层内与层间效应

4.1 为什么中心化会改变研究问题

对群组内个体变量 XijX_{ij}：

Xij=(Xij−X‾j)+X‾j, X_{ij}=(X_{ij}-\bar X_j)+\bar X_j, Yij=β0+βW(Xij−X‾j)+βBX‾j+u0j+εij. Y_{ij}=\beta_0+\beta_W(X_{ij}-\bar X_j)+ \beta_B\bar X_j+u_{0j}+\varepsilon_{ij}.

  • βW\beta_W：同一群组内两名个体相差一单位时的关系；
  • βB\beta_B：群组平均水平相差一单位时的关系；
  • βB−βW\beta_B-\beta_W：常称 contextual contrast，但仍需额外假设才能作因果解释。

总体均值中心化主要改变零点解释和数值稳定性，不会自动分离层内与层间关系。群组均值中心化有助于分解问题，但也不会自动消除所有群组层混杂。

一个系数可能混合两个方向不同的问题 把时间变化压力直接放入模型，会混合“同一个人今天比平时压力更高”和“这个人长期比别人压力更高”。同理，学生家庭 SES 与学校整体 SES 也不应由一个未分解系数概括。

5 社会学应用：学生嵌套于学校

5.1 模拟 SES、学校支持与成绩

set.seed(20260822)
J <- 70L
n_j <- sample(24:36, J, replace = TRUE)
school <- sprintf("S%02d", seq_len(J))
school_ses <- rnorm(J, 0, 0.8)
school_ses_gc <- school_ses - mean(school_ses)
support_gc <- rnorm(J, 0, 1)
support_gc <- support_gc - mean(support_gc)
rho <- 0.25
z0 <- rnorm(J)
z1 <- rnorm(J)
b0 <- 5.0 * z0
b1 <- 1.5 * (rho * z0 + sqrt(1 - rho^2) * z1)

soc <- do.call(rbind, lapply(seq_len(J), function(j) {
  n <- n_j[j]
  ses_wc <- rnorm(n, 0, 0.9)
  ses_wc <- ses_wc - mean(ses_wc)
  ses <- school_ses[j] + ses_wc
  achievement <- 70 +
    2.8 * ses_wc + 5.0 * school_ses_gc[j] +
    2.0 * support_gc[j] - 1.2 * ses_wc * support_gc[j] +
    b0[j] + b1[j] * ses_wc + rnorm(n, 0, 7)
  data.frame(
    school = school[j],
    ses = ses,
    school_support_gc = support_gc[j],
    achievement = achievement
  )
}))
soc$school_ses <- ave(soc$ses, soc$school, FUN = mean)
soc$ses_wc <- soc$ses - soc$school_ses
soc$school_ses_gc <- soc$school_ses - mean(soc$school_ses)

soc_audit <- data.frame(
  学生数 = nrow(soc),
  学校数 = length(unique(soc$school)),
  最小学校人数 = min(table(soc$school)),
  最大学校人数 = max(table(soc$school)),
  `ses_wc 校内均值最大绝对值` = max(abs(tapply(soc$ses_wc, soc$school, mean))),
  check.names = FALSE
)

knitr::kable(soc_audit, digits = 5, caption = "学生—学校模拟数据与中心化审计")
学生—学校模拟数据与中心化审计
学生数 学校数 最小学校人数 最大学校人数 ses_wc 校内均值最大绝对值
2098 70 24 36 0

ses_wc 在每所学校内均值为零，回答个体相对本校同伴的位置；这里的 school_ses_gc 是学校平均 SES 相对按学生人数加权的总均值的位置。若研究问题要求每所学校等权，应先提取每校唯一均值再中心化，并在分析计划中说明。学校支持则在学校层按学校等权进行总体均值中心化。

school_means <- aggregate(
  cbind(achievement, school_ses) ~ school,
  data = soc,
  FUN = mean
)
old_par <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 3, 1))
plot(
  soc$ses_wc,
  soc$achievement,
  pch = 16,
  cex = 0.45,
  col = grDevices::adjustcolor(palette_ml["teal"], alpha.f = 0.30),
  xlab = "学生 SES − 本校平均 SES",
  ylab = "成绩",
  main = "校内关系",
  las = 1
)
abline(stats::lm(achievement ~ ses_wc, data = soc), col = palette_ml["navy"], lwd = 2)
plot(
  school_means$school_ses,
  school_means$achievement,
  pch = 17,
  col = palette_ml["orange"],
  xlab = "学校平均 SES",
  ylab = "学校平均成绩",
  main = "校间关系",
  las = 1
)
abline(
  stats::lm(achievement ~ school_ses, data = school_means),
  col = palette_ml["vermillion"],
  lwd = 2
)
两个散点图。左侧有两千多个学生点，显示校内 SES 与成绩的正相关；右侧七十个学校均值点显示更陡的学校间正相关。

学生 SES 的校内与校间关系。左图展示学生相对本校平均 SES 与成绩，右图展示学校平均 SES 与学校平均成绩；两者是不同层级的问题。 

par(old_par)

5.2 随机斜率与跨层交互模型

随机截距模型假定 SES 的校内斜率在每所学校相同。允许斜率变化时：

Yij=β0+β1Xij+u0j+u1jXij+εij, Y_{ij}=\beta_0+\beta_1X_{ij}+u_{0j}+u_{1j}X_{ij}+\varepsilon_{ij}, (u0ju1j)∼N[(00),(τ00τ01τ01τ11)]. \begin{pmatrix}u_{0j}\\u_{1j}\end{pmatrix} \sim N\left[ \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}\tau_{00}&\tau_{01}\\\tau_{01}&\tau_{11}\end{pmatrix} \right].

β1\beta_1 是平均校内斜率，τ11\tau_{11} 描述学校间斜率异质性，τ01\tau_{01} 描述截距与斜率如何共同变化。加入学校支持 ZjZ_j 的跨层交互后，平均 SES 斜率为 β1+β3Zj\beta_1+\beta_3Z_j。

m_soc_empty <- lme4::lmer(
  achievement ~ 1 + (1 | school),
  data = soc,
  REML = TRUE
)

m_soc <- lme4::lmer(
  achievement ~ ses_wc + school_ses_gc + school_support_gc +
    ses_wc:school_support_gc + (1 + ses_wc | school),
  data = soc,
  REML = TRUE,
  control = lme4::lmerControl(optimizer = "bobyqa")
)

soc_labels <- c(
  `(Intercept)` = "平均 SES、平均支持学校的截距",
  ses_wc = "学生 SES 的校内效应",
  school_ses_gc = "学校平均 SES 的校间效应",
  school_support_gc = "学校支持的主效应",
  `ses_wc:school_support_gc` = "校内 SES × 学校支持"
)
soc_fixed <- lmm_fixed_table(m_soc, soc_labels)
soc_empty_icc <- icc_random_intercept(m_soc_empty, "school")

soc_vc_raw <- as.data.frame(lme4::VarCorr(m_soc))
soc_vc_table <- data.frame(
  群组 = c("学校", "学校", "学校", "残差"),
  成分 = c("随机截距", "SES 随机斜率", "截距–斜率相关", "学生层残差"),
  数值 = c(
    soc_vc_raw$sdcor[
      soc_vc_raw$grp == "school" & soc_vc_raw$var1 == "(Intercept)" &
        is.na(soc_vc_raw$var2)
    ],
    soc_vc_raw$sdcor[
      soc_vc_raw$grp == "school" & soc_vc_raw$var1 == "ses_wc" &
        is.na(soc_vc_raw$var2)
    ],
    soc_vc_raw$sdcor[
      soc_vc_raw$grp == "school" & !is.na(soc_vc_raw$var2)
    ],
    stats::sigma(m_soc)
  ),
  check.names = FALSE
)

knitr::kable(
  soc_fixed,
  digits = 3,
  caption = "学生—学校随机斜率模型的固定效应"
)
学生—学校随机斜率模型的固定效应
变量 估计值 标准误 95% CI 下限 95% CI 上限 t 值
(Intercept) 平均 SES、平均支持学校的截距 70.88 0.644 69.62 72.143 110.06
ses_wc 学生 SES 的校内效应 2.77 0.260 2.26 3.274 10.64
school_ses_gc 学校平均 SES 的校间效应 4.32 0.719 2.91 5.732 6.01
school_support_gc 学校支持的主效应 2.78 0.712 1.38 4.173 3.90
ses_wc:school_support_gc 校内 SES × 学校支持 -1.01 0.285 -1.57 -0.454 -3.55 
knitr::kable(
  soc_vc_table,
  digits = 3,
  caption = "学生—学校模型的随机效应标准差、相关与残差"
)
学生—学校模型的随机效应标准差、相关与残差
群组 成分 数值
学校 随机截距 5.233
学校 SES 随机斜率 1.598
学校 截距–斜率相关 0.287
残差 学生层残差 6.994

校内 SES 关系估计为每增加 1 单位，成绩平均增加 2.77 分；学校平均 SES 的校间关系为 4.32 分。两者不同，不能用一个原始 SES 系数代替。跨层交互为 -1.01：学校支持越高，SES 的正向校内斜率在本模拟中越弱。它是效应修饰描述，不证明提高支持必然改变 SES 机制。

空模型学校 ICC 为 0.437。随机斜率标准差为 1.598，说明校内 SES 关系存在学校间异质性；截距–斜率相关为 0.287，其解释依赖 ses_wc=0 的零点。

soc_re <- lme4::ranef(m_soc)$school
selected_schools <- rownames(soc_re)[round(seq(1, nrow(soc_re), length.out = 12))]
school_profile <- unique(
  soc[c("school", "school_ses_gc", "school_support_gc")]
)
ses_grid <- seq(-2, 2, length.out = 80)
soc_beta <- lme4::fixef(m_soc)
school_predictions <- stats::setNames(
  lapply(selected_schools, function(s) {
    profile_s <- school_profile[school_profile$school == s, ]
    support_s <- profile_s$school_support_gc
    school_ses_s <- profile_s$school_ses_gc
    soc_beta["(Intercept)"] +
      soc_beta["school_ses_gc"] * school_ses_s +
      soc_beta["school_support_gc"] * support_s +
      soc_re[s, "(Intercept)"] +
      (soc_beta["ses_wc"] +
         soc_beta["ses_wc:school_support_gc"] * support_s +
         soc_re[s, "ses_wc"]) * ses_grid
  }),
  selected_schools
)
average_school_prediction <-
  soc_beta["(Intercept)"] + soc_beta["ses_wc"] * ses_grid
left_ylim <- grDevices::extendrange(
  c(unlist(school_predictions), average_school_prediction),
  f = 0.05
)
old_par <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 3, 1))

plot(
  NA,
  xlim = range(ses_grid),
  ylim = left_ylim,
  xlab = "学生 SES − 本校平均 SES",
  ylab = "条件预测成绩",
  main = "学校特异随机斜率",
  las = 1
)
for (s in selected_schools) {
  lines(
    ses_grid,
    school_predictions[[s]],
    col = grDevices::adjustcolor(palette_ml["blue"], alpha.f = 0.45),
    lwd = 1.2
  )
}
lines(
  ses_grid,
  average_school_prediction,
  col = palette_ml["vermillion"],
  lwd = 2.6
)

support_sd <- stats::sd(
  soc$school_support_gc[!duplicated(soc$school)]
)
support_values <- c(-support_sd, 0, support_sd)
support_colors <- c(
  palette_ml["orange"], palette_ml["navy"], palette_ml["teal"]
)
support_predictions <- vapply(
  support_values,
  function(z) {
    soc_beta["(Intercept)"] + soc_beta["school_support_gc"] * z +
      (soc_beta["ses_wc"] +
         soc_beta["ses_wc:school_support_gc"] * z) * ses_grid
  },
  numeric(length(ses_grid))
)
right_ylim <- grDevices::extendrange(support_predictions, f = 0.05)
plot(
  NA,
  xlim = range(ses_grid),
  ylim = right_ylim,
  xlab = "学生 SES − 本校平均 SES",
  ylab = "固定部分预测成绩",
  main = "学校支持的跨层修饰",
  las = 1
)
for (k in seq_along(support_values)) {
  lines(
    ses_grid,
    support_predictions[, k],
    col = support_colors[k],
    lwd = 2.4
  )
}
legend(
  "topleft",
  legend = c("支持 = −1 SD", "支持 = 平均", "支持 = +1 SD"),
  col = support_colors,
  lwd = 2.4,
  bty = "n"
)
两个面板。左侧多条学校特异成绩线具有不同截距和斜率；右侧三条固定预测线显示学校支持越高，SES 斜率越平缓。

学生—学校随机斜率与跨层交互。左图展示部分学校的条件 SES 斜率，右图展示学校支持为低、平均和高水平时的固定部分预测关系。 

par(old_par)

5.3 随机斜率下相关性随协变量改变

当两名同校学生的中心化 SES 分别为 x1,x2x_1,x_2 时，随机效应诱导的协方差为：

Cov⁡(Y1,Y2)=τ00+(x1+x2)τ01+x1x2τ11. \operatorname{Cov}(Y_1,Y_2)= \tau_{00}+(x_1+x_2)\tau_{01}+x_1x_2\tau_{11}.

因此不再有一个对所有 SES 都相同的 ICC。

G_soc <- as.matrix(lme4::VarCorr(m_soc)$school)
sigma2_soc <- stats::sigma(m_soc)^2

conditional_corr <- function(x1, x2, G, sigma2) {
  z1 <- c(1, x1)
  z2 <- c(1, x2)
  covariance <- drop(t(z1) %*% G %*% z2)
  variance1 <- drop(t(z1) %*% G %*% z1) + sigma2
  variance2 <- drop(t(z2) %*% G %*% z2) + sigma2
  covariance / sqrt(variance1 * variance2)
}

soc_corr_table <- data.frame(
  学生一SES = c(0, -1),
  学生二SES = c(0, 1),
  同校模型相关 = c(
    conditional_corr(0, 0, G_soc, sigma2_soc),
    conditional_corr(-1, 1, G_soc, sigma2_soc)
  ),
  check.names = FALSE
)

knitr::kable(
  soc_corr_table,
  digits = 3,
  caption = "随机斜率模型在不同 SES 组合下的同校相关"
)
随机斜率模型在不同 SES 组合下的同校相关
学生一SES 学生二SES 同校模型相关
0 0 0.359
-1 1 0.315
检验你的理解：为什么随机截距不能自动替代随机斜率？ 随机截距只允许学校平均水平不同，仍强制 SES 斜率完全相同。如果科学问题和重复信息支持斜率异质性，只放随机截距会错误限制协方差结构，并可能使固定斜率的不确定性过小。

6 心理学应用：时点—个体—治疗师

6.1 三层纵向数据

160 名参与者分别由 20 名治疗师服务，每人记录 7 个治疗 session。时间从 0 开始，使截距表示基线症状。压力被分成个人内偏离 stress_wp 与个人平均压力 stress_pm_gc。

set.seed(20260823)
K <- 20L
therapists <- sprintf("T%02d", seq_len(K))
therapist_u <- rnorm(K, 0, 1.2)
persons_per_therapist <- 8L
session0 <- 0:6
psy_list <- list()
person_index <- 0L

for (k in seq_len(K)) {
  intervention <- sample(rep(0:1, each = persons_per_therapist / 2))
  for (p in seq_len(persons_per_therapist)) {
    person_index <- person_index + 1L
    z0p <- rnorm(1)
    z1p <- rnorm(1)
    p_u0 <- 3.0 * z0p
    p_u1 <- 0.45 * (-0.20 * z0p + sqrt(1 - 0.20^2) * z1p)
    stress_pm <- rnorm(1, 0, 0.9)
    stress_wp <- rnorm(length(session0), 0, 0.8)
    stress_wp <- stress_wp - mean(stress_wp)
    symptoms <- 22 - 0.90 * session0 -
      0.20 * intervention[p] - 0.35 * session0 * intervention[p] +
      1.20 * stress_wp + 2.00 * stress_pm +
      therapist_u[k] + p_u0 + p_u1 * session0 +
      rnorm(length(session0), 0, 2.5)
    psy_list[[person_index]] <- data.frame(
      therapist = therapists[k],
      person = sprintf("P%03d", person_index),
      session0 = session0,
      intervention = intervention[p],
      stress = stress_pm + stress_wp,
      symptoms = symptoms
    )
  }
}
psy <- do.call(rbind, psy_list)
psy$stress_pm <- ave(psy$stress, psy$person, FUN = mean)
psy$stress_wp <- psy$stress - psy$stress_pm
psy$stress_pm_gc <- psy$stress_pm - mean(psy$stress_pm)
psy$intervention <- factor(
  psy$intervention,
  0:1,
  c("Control", "Intervention")
)

psy_audit <- data.frame(
  观测数 = nrow(psy),
  个体数 = length(unique(psy$person)),
  治疗师数 = length(unique(psy$therapist)),
  每人时点数 = paste(range(table(psy$person)), collapse = "–"),
  每名治疗师人数 = paste(range(table(unique(psy[c("person", "therapist")])$therapist)), collapse = "–"),
  check.names = FALSE
)

knitr::kable(psy_audit, caption = "心理学三层纵向数据审计")
心理学三层纵向数据审计
观测数 个体数 治疗师数 每人时点数 每名治疗师人数
1120 160 20 7–7 8–8

6.2 三层随机截距与个体时间随机斜率

Ytij=β0+β1Timetij+β2Aij+β3TimetijAij+βWStresswithin+βBStressbetween+v0j+u0ij+u1ijTimetij+εtij. Y_{tij}=\beta_0+\beta_1Time_{tij}+\beta_2A_{ij}+ \beta_3Time_{tij}A_{ij}+\beta_WStress_{within}+ \beta_BStress_{between}+v_{0j}+u_{0ij}+u_{1ij}Time_{tij}+\varepsilon_{tij}. 

m_psy <- lme4::lmer(
  symptoms ~ session0 * intervention + stress_wp + stress_pm_gc +
    (1 + session0 | person) + (1 | therapist),
  data = psy,
  REML = TRUE,
  control = lme4::lmerControl(optimizer = "bobyqa")
)

psy_labels <- c(
  `(Intercept)` = "对照组基线症状",
  session0 = "对照组每 session 平均变化",
  interventionIntervention = "干预组基线差",
  stress_wp = "压力的个人内效应",
  stress_pm_gc = "个人平均压力的个人间效应",
  `session0:interventionIntervention` = "干预组额外 session 变化"
)
psy_fixed <- lmm_fixed_table(m_psy, psy_labels)

psy_vc_raw <- as.data.frame(lme4::VarCorr(m_psy))
psy_vc_table <- data.frame(
  成分 = c("个体随机截距 SD", "个体时间随机斜率 SD", "截距–斜率相关", "治疗师截距 SD", "残差 SD"),
  数值 = c(
    psy_vc_raw$sdcor[
      psy_vc_raw$grp == "person" & psy_vc_raw$var1 == "(Intercept)" &
        is.na(psy_vc_raw$var2)
    ],
    psy_vc_raw$sdcor[
      psy_vc_raw$grp == "person" & psy_vc_raw$var1 == "session0" &
        is.na(psy_vc_raw$var2)
    ],
    psy_vc_raw$sdcor[
      psy_vc_raw$grp == "person" & !is.na(psy_vc_raw$var2)
    ],
    psy_vc_raw$sdcor[
      psy_vc_raw$grp == "therapist" & is.na(psy_vc_raw$var2)
    ],
    stats::sigma(m_psy)
  ),
  check.names = FALSE
)

knitr::kable(
  psy_fixed,
  digits = 3,
  caption = "时点—个体—治疗师三层模型固定效应"
)
时点—个体—治疗师三层模型固定效应
变量 估计值 标准误 95% CI 下限 95% CI 上限 t 值
(Intercept) 对照组基线症状 22.335 0.524 21.307 23.363 42.59
session0 对照组每 session 平均变化 -0.857 0.062 -0.979 -0.735 -13.74
interventionIntervention 干预组基线差 -0.621 0.531 -1.661 0.419 -1.17
stress_wp 压力的个人内效应 1.326 0.104 1.123 1.529 12.78
stress_pm_gc 个人平均压力的个人间效应 1.654 0.289 1.088 2.220 5.73
session0:interventionIntervention 干预组额外 session 变化 -0.403 0.088 -0.576 -0.230 -4.57 
knitr::kable(
  psy_vc_table,
  digits = 3,
  caption = "心理学三层模型的随机效应与残差"
)
心理学三层模型的随机效应与残差
成分 数值
个体随机截距 SD 2.869
个体时间随机斜率 SD 0.295
截距–斜率相关 0.043
治疗师截距 SD 1.642
残差 SD 2.508

对照组每个 session 症状平均改变 -0.86 分；干预组在此基础上每个 session 额外改变 -0.4 分。当天压力比个人通常水平高 1 单位时，症状平均高 1.33 分；长期平均压力高 1 单位的人之间相差 1.65 分。这两个系数回答不同问题。

只有 20 名治疗师，治疗师方差和任何治疗师排名都具有较大不确定性。随机治疗师截距处理依赖结构，但不会证明治疗师差异由治疗质量造成。

selected_people <- unique(psy$person)[round(seq(1, length(unique(psy$person)), length.out = 12))]
psy_selected <- psy[psy$person %in% selected_people, ]
psy_beta <- lme4::fixef(m_psy)

plot(
  NA,
  xlim = range(psy$session0),
  ylim = range(psy_selected$symptoms),
  xlab = "Session（0 = 基线）",
  ylab = "症状评分",
  main = "个体轨迹与固定部分平均轨迹",
  las = 1
)
for (id in selected_people) {
  d <- psy_selected[psy_selected$person == id, ]
  color <- if (d$intervention[1] == "Intervention") {
    grDevices::adjustcolor(palette_ml["teal"], alpha.f = 0.45)
  } else {
    grDevices::adjustcolor(palette_ml["orange"], alpha.f = 0.45)
  }
  lines(d$session0, d$symptoms, col = color, lwd = 1.2)
  points(d$session0, d$symptoms, col = color, pch = 16, cex = 0.45)
}
session_grid <- 0:6
fixed_control <- psy_beta["(Intercept)"] + psy_beta["session0"] * session_grid
fixed_intervention <- psy_beta["(Intercept)"] +
  psy_beta["interventionIntervention"] +
  (psy_beta["session0"] +
     psy_beta["session0:interventionIntervention"]) * session_grid
lines(session_grid, fixed_control, col = palette_ml["orange"], lwd = 3)
lines(session_grid, fixed_intervention, col = palette_ml["teal"], lwd = 3)
legend(
  "topright",
  legend = c("对照固定轨迹", "干预固定轨迹"),
  col = c(palette_ml["orange"], palette_ml["teal"]),
  lwd = 3,
  bty = "n"
)
十二条个体症状折线随七次 session 总体下降但斜率不同。两条粗线显示干预组平均下降比对照组更快。

部分参与者的症状轨迹及模型固定部分的平均轨迹。细线为观察到的个体轨迹，粗线在压力取平均值时比较对照与干预的固定部分变化。

随机时间斜率允许个体变化速度不同，但不自动消除相邻 session 残差的 AR(1) 相关。若残差仍有明显时间序列结构，应使用支持相应相关结构的方法并做敏感性分析。

7 医学应用二：患者嵌套于医院的二分类结局

7.1 为什么需要 logistic GLMM

二分类感染结局满足：

Yij∣uj∼Bernoulli(pij),logit⁡(pij)=XijTβ+uj. Y_{ij}\mid u_j\sim Bernoulli(p_{ij}),\qquad \operatorname{logit}(p_{ij})=X_{ij}^{T}\beta+u_j.

这里 exp⁡(β)\exp(\beta) 是给定医院随机效应与模型协变量后的 cluster-specific 条件 OR，通常不等于对医院分布平均后的 marginal OR。

set.seed(20260829)
J <- 56L
n_j <- sample(36:54, J, replace = TRUE)
hospital <- sprintf("H%02d", seq_len(J))
checklist_j <- sample(rep(0:1, each = J / 2))
h_u <- rnorm(J, 0, 0.65)

bin <- do.call(rbind, lapply(seq_len(J), function(j) {
  n <- n_j[j]
  age <- pmin(pmax(rnorm(n, 62, 12), 18), 90)
  emergency <- rbinom(n, 1, 0.28)
  eta <- -2.25 - 0.55 * checklist_j[j] +
    0.70 * emergency + 0.22 * ((age - 60) / 10) + h_u[j]
  data.frame(
    hospital = hospital[j],
    checklist = checklist_j[j],
    emergency = emergency,
    age_c10 = (age - 60) / 10,
    infection = rbinom(n, 1, plogis(eta))
  )
}))
bin$checklist <- factor(bin$checklist, 0:1, c("Usual care", "Checklist"))
bin$emergency <- factor(bin$emergency, 0:1, c("Elective", "Emergency"))

bin_audit <- data.frame(
  患者数 = nrow(bin),
  医院数 = length(unique(bin$hospital)),
  感染事件 = sum(bin$infection),
  感染比例 = pct(mean(bin$infection)),
  医院人数范围 = paste(range(table(bin$hospital)), collapse = "–"),
  check.names = FALSE
)

knitr::kable(bin_audit, caption = "医院感染二分类模拟数据审计")
医院感染二分类模拟数据审计
患者数 医院数 感染事件 感染比例 医院人数范围
2501 56 298 11.9% 36–53 
hospital_rates <- aggregate(infection ~ hospital, data = bin, FUN = mean)
hospital_rates$n <- as.numeric(table(bin$hospital)[hospital_rates$hospital])
hospital_rates <- hospital_rates[order(hospital_rates$infection), ]

plot(
  hospital_rates$infection,
  seq_len(nrow(hospital_rates)),
  pch = 16,
  cex = 0.7 + 1.2 * sqrt(hospital_rates$n / max(hospital_rates$n)),
  col = palette_ml["blue"],
  yaxt = "n",
  xlab = "观察感染比例",
  ylab = "按感染比例排序的医院",
  main = "医院原始感染比例",
  las = 1
)
abline(v = mean(bin$infection), lty = 2, col = palette_ml["vermillion"])
五十六家医院的感染比例点按数值排序，围绕总体感染率虚线分布；部分小样本医院位于两端。

各医院观察到的感染比例及总体比例。医院样本量和病例结构不同，原始比例不应直接解释为医院质量排名。

7.2 条件 OR、latent ICC 与 median odds ratio

m_bin <- lme4::glmer(
  infection ~ checklist + emergency + age_c10 + (1 | hospital),
  data = bin,
  family = stats::binomial,
  nAGQ = 1,
  control = lme4::glmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 2e5)
  )
)

bin_vc <- as.data.frame(lme4::VarCorr(m_bin))
tau2_bin <- bin_vc$vcov[
  bin_vc$grp == "hospital" &
    bin_vc$var1 == "(Intercept)" &
    is.na(bin_vc$var2)
]
stopifnot(length(tau2_bin) == 1L)
icc_latent <- tau2_bin / (tau2_bin + pi^2 / 3)
mor <- exp(stats::qnorm(0.75) * sqrt(2 * tau2_bin))
pearson <- residuals(m_bin, type = "pearson")
pearson_ratio <- sum(pearson^2) /
  (stats::nobs(m_bin) - length(lme4::fixef(m_bin)))

bin_coef <- coef(summary(m_bin))
bin_or <- wald_or(m_bin)
bin_labels <- c(
  `(Intercept)` = "常规护理、择期、60 岁",
  checklistChecklist = "感染检查清单 vs 常规护理",
  emergencyEmergency = "急诊 vs 择期",
  age_c10 = "年龄每增加 10 岁"
)
bin_results <- data.frame(
  变量 = unname(bin_labels[rownames(bin_coef)]),
  条件OR = bin_or[, "OR"],
  `95% CI 下限` = bin_or[, "lower"],
  `95% CI 上限` = bin_or[, "upper"],
  `p 值` = vapply(bin_coef[, "Pr(>|z|)"], format_p, character(1)),
  check.names = FALSE
)

bin_heterogeneity <- data.frame(
  指标 = c("医院随机截距方差", "latent-scale ICC", "median odds ratio", "Pearson 比率"),
  数值 = c(tau2_bin, icc_latent, mor, pearson_ratio),
  check.names = FALSE
)

knitr::kable(
  bin_results,
  digits = 3,
  caption = "医院感染 logistic GLMM 的条件 OR 与 Wald 区间"
)
医院感染 logistic GLMM 的条件 OR 与 Wald 区间
变量 条件OR 95% CI 下限 95% CI 上限 p 值
(Intercept) 常规护理、择期、60 岁 0.126 0.094 0.169 <0.001
checklistChecklist 感染检查清单 vs 常规护理 0.571 0.380 0.856 0.00676
emergencyEmergency 急诊 vs 择期 1.562 1.197 2.037 0.001
age_c10 年龄每增加 10 岁 1.254 1.130 1.393 <0.001 
knitr::kable(
  bin_heterogeneity,
  digits = 3,
  caption = "医院异质性与二分类模型诊断摘要"
)
医院异质性与二分类模型诊断摘要
指标 数值
医院随机截距方差 0.347
latent-scale ICC 0.095
median odds ratio 1.754
Pearson 比率 0.905

模型显式使用 nAGQ = 1，即 lme4::glmer() 的 Laplace 近似。正式分析应报告积分近似与 optimizer；在事件稀少、群组很小或随机效应较大时，可增加自适应 Gauss–Hermite quadrature 节点并检查结论敏感性。

检查清单的条件 OR 为 0.571（95% CI：0.38 至 0.856）。在相同年龄、急诊状态和相同潜在医院倾向下，检查清单组的感染 odds 较低。

简单 logistic 随机截距模型的 latent ICC 为：

ICClatent≈τ00τ00+π2/3. ICC_{latent}\approx\frac{\tau_{00}}{\tau_{00}+\pi^2/3}.

本例为 0.095，它属于潜在 logistic 尺度，不是观察感染 0/1 方差中可直接解释的百分比。Median odds ratio 为 1.754：随机选择两个风险不同的医院，把相同患者从低风险医院移到高风险医院时，医院效应的中位 odds 倍数约为此值。

Pearson 比率为 0.905，没有显示明显过度离散，但单一比率不能排除函数形式错误、零膨胀或遗漏层级。

7.3 将医院效应设为零不等于对医院异质性积分

beta <- lme4::fixef(m_bin)
eta_usual <- unname(beta["(Intercept)"])
eta_check <- eta_usual + unname(beta["checklistChecklist"])
p_typical <- stats::plogis(c(usual = eta_usual, checklist = eta_check))

integrand <- function(u, eta, sd) {
  stats::plogis(eta + u) * stats::dnorm(u, 0, sd)
}

p_marginal <- c(
  usual = stats::integrate(
    integrand,
    -Inf,
    Inf,
    eta = eta_usual,
    sd = sqrt(tau2_bin)
  )$value,
  checklist = stats::integrate(
    integrand,
    -Inf,
    Inf,
    eta = eta_check,
    sd = sqrt(tau2_bin)
  )$value
)

probability_table <- data.frame(
  估计尺度 = c(
    "60 岁择期参考患者、典型医院：u=0",
    "60 岁择期参考患者、对医院随机效应积分"
  ),
  常规护理风险 = c(p_typical["usual"], p_marginal["usual"]),
  检查清单风险 = c(p_typical["checklist"], p_marginal["checklist"]),
  风险差 = c(diff(p_typical), diff(p_marginal)),
  check.names = FALSE
)

knitr::kable(
  probability_table,
  digits = 3,
  caption = "参考患者在典型医院与对医院随机效应积分后的感染概率"
)
参考患者在典型医院与对医院随机效应积分后的感染概率
估计尺度 常规护理风险 检查清单风险 风险差
60 岁择期参考患者、典型医院：u=0 0.112 0.067 -0.045
60 岁择期参考患者、对医院随机效应积分 0.125 0.077 -0.048

predict(m_bin, re.form = NA, type = "response") 把医院随机效应设为零，得到典型医院的 fixed-part 概率；它不会自动对医院分布积分。上表第二行仅对医院随机截距积分，仍以 60 岁、择期入院者为参考。若目标是研究人群的总体边际风险，还需按目标人群的年龄、急诊状态等协变量分布进一步标准化。

prob_matrix <- rbind(
  `参考患者：典型医院 u=0` = p_typical,
  `参考患者：医院效应积分` = p_marginal
)
barplot(
  t(prob_matrix),
  beside = TRUE,
  col = c(palette_ml["orange"], palette_ml["teal"]),
  ylim = c(0, max(prob_matrix) * 1.25),
  ylab = "预测感染概率",
  main = "条件固定部分与边际化概率",
  las = 1
)
legend(
  "topright",
  legend = c("常规护理", "检查清单"),
  fill = c(palette_ml["orange"], palette_ml["teal"]),
  bty = "n"
)
分组柱状图比较常规护理和检查清单。两种方法均显示检查清单风险较低；对医院随机效应积分后的两根柱均略高于典型医院概率。

60 岁择期参考患者在典型医院与对医院随机效应积分后的感染概率。非线性 logit 链接使积分结果不同于把随机效应设为零的概率。

检验你的理解：条件 OR=0.57 能否写成感染风险降低 43%？ 不能。OR 比较 odds，不是风险比；而且这里是给定医院随机效应的条件 OR。应另行计算目标人群和协变量分布下的标准化风险、风险比或风险差，并明确是否对随机效应积分。

8 更复杂的群组结构

8.1 嵌套、交叉分类与多重成员关系不是同一语法

若学校严格嵌套于社区，可写作 (1 | community/school)；若学生的居住社区与学校互相交叉，应分别估计两套随机截距。把交叉分类强行写成嵌套模型，会把方差归给错误的社会结构。

# 学校严格嵌套于社区；school_id 只需在社区内唯一。
m_nested <- lme4::lmer(
  score ~ student_ses + school_resources +
    (1 | community_id/school_id),
  data = education_data
)

# 学校和居住社区交叉分类。
m_crossed <- lme4::lmer(
  score ~ student_ses + school_resources + neighborhood_deprivation +
    (1 | school_id) + (1 | neighborhood_id),
  data = education_data
)

# 多重成员关系：若一名患者由多名治疗师共同服务，需要明确每名
# 治疗师的成员权重；普通单一 (1 | therapist) 不能表达这种结构。

多重成员模型通常需要每个上层成员的权重及支持该结构的软件。若患者更换治疗师，只保留“最后一名治疗师”会丢失真实依赖过程。

9 ML、REML 与推断

9.1 估计方法服务于比较目标

fitting_reference <- data.frame(
  情形 = c(
    "LMM 最终方差成分与固定效应",
    "比较不同固定效应的 LMM",
    "相同固定效应、比较协方差结构",
    "logistic／Poisson GLMM",
    "少量高层群组"
  ),
  建议 = c(
    "常用 REML，并报告固定与随机部分",
    "用 ML 重拟合后再比较",
    "可用 REML，但边界检验需谨慎",
    "使用 maximum likelihood；说明积分近似与 optimizer",
    "优先效应与区间；考虑 bootstrap 或小样本校正"
  ),
  原因 = c(
    "REML 减少方差成分有限样本偏差",
    "不同固定效应模型的 REML likelihood 不可直接比较",
    "随机方差为零处在参数空间边界",
    "非高斯模型没有普通 LMM 的 REML",
    "渐近 z/t 近似可能过于乐观"
  ),
  check.names = FALSE
)

knitr::kable(fitting_reference, caption = "ML、REML 与不确定性方法的选择")
ML、REML 与不确定性方法的选择
情形 建议 原因
LMM 最终方差成分与固定效应 常用 REML，并报告固定与随机部分 REML 减少方差成分有限样本偏差
比较不同固定效应的 LMM 用 ML 重拟合后再比较 不同固定效应模型的 REML likelihood 不可直接比较
相同固定效应、比较协方差结构 可用 REML，但边界检验需谨慎 随机方差为零处在参数空间边界
logistic／Poisson GLMM 使用 maximum likelihood；说明积分近似与 optimizer 非高斯模型没有普通 LMM 的 REML
少量高层群组 优先效应与区间；考虑 bootstrap 或小样本校正 渐近 z/t 近似可能过于乐观

lme4::lmer() 默认不提供固定效应 p 值，因为有限样本自由度并不唯一。教程使用估计值、Wald 95% CI 和 t 值描述结果；正式研究可预先指定 Satterthwaite、Kenward–Roger、profile likelihood 或 parametric bootstrap。不要在看到结果后选择最有利的区间方法。

随机效应方差的零假设位于边界，普通似然比检验的卡方近似可能不准确。随机结构应首先来自数据层级、重复测量安排和科学问题，而不是逐步寻找最小 p 值。

更多低层观测不能替代更多高层群组 每家医院增加患者能改善该医院均值，但医院层政策、治疗或资源效应的精度主要受医院数限制。报告总行数时必须同时报告每一层的单位数和群组大小分布。

9.2 条件预测、典型群组与新群组

# 已观察个体／群组的条件预测：包含其估计随机效应。
predict(m_psy, newdata = existing_people, re.form = NULL)

# 把所有随机效应设为零：固定部分／典型群组预测。
predict(m_psy, newdata = prediction_grid, re.form = NA)

# 新个体由既有治疗师服务：保留治疗师效应，不包含未知的个体效应。
predict(
  m_psy,
  newdata = new_people_existing_therapist,
  re.form = ~(1 | therapist),
  allow.new.levels = TRUE
)

# 新个体和新治疗师：两层随机效应都未知，使用 fixed part。
predict(
  m_psy,
  newdata = new_people_new_therapist,
  re.form = NA,
  allow.new.levels = TRUE
)

re.form = NA 会同时清除模型中的所有随机效应；它不能表示“新个体、既有治疗师”这种部分已知结构。后者需要像示例一样只保留治疗师项。LMM identity link 下，对零均值随机效应平均后的均值与 fixed part 对齐；GLMM 的非线性链接下，把随机效应设为零与对其分布积分通常不同。预测还应区分均值置信区间、已有个体的条件预测区间和新群组预测区间。

10 诊断：模型能运行不等于模型可信

10.1 统一检查四个主模型

diagnostic_flags <- data.frame(
  模型 = c(
    "医学连续 LMM", "社会学随机斜率 LMM",
    "心理学三层 LMM", "医学感染 GLMM"
  ),
  奇异拟合 = c(
    lme4::isSingular(m_med),
    lme4::isSingular(m_soc),
    lme4::isSingular(m_psy),
    lme4::isSingular(m_bin)
  ),
  优化器报告收敛 = vapply(
    list(m_med, m_soc, m_psy, m_bin),
    function(model) {
      optimizer_code <- model@optinfo$conv$opt
      all(optimizer_code == 0) &&
        is.null(model@optinfo$conv$lme4$messages)
    },
    logical(1)
  ),
  check.names = FALSE
)

knitr::kable(diagnostic_flags, caption = "四个教学模型的奇异性与收敛标志")
四个教学模型的奇异性与收敛标志
模型 奇异拟合 优化器报告收敛
医学连续 LMM FALSE TRUE
社会学随机斜率 LMM FALSE TRUE
心理学三层 LMM FALSE TRUE
医学感染 GLMM FALSE TRUE

四个教学模型均成功收敛且不是 singular fit，但这只是最低要求。完整诊断还包括：

  • 每层 ID、重复记录和群组大小是否正确；
  • 连续 predictor 在拟合 random slope 的群组内是否有足够变化；
  • 固定部分是否存在非线性、遗漏交互或异方差；
  • observation 与 cluster 两个层级的影响点；
  • random-effect 分布是否有严重偏离；
  • 二分类模型是否过度离散、零膨胀或遗漏层级；
  • 纵向残差是否仍存在时间序列相关；
  • 结论是否由少数极大或极小群组驱动。

10.1.1 收敛 warning 与 singular fit

收敛 warning 可能来自尺度差异、梯度、过度复杂结构或弱信息。建议先核对数据、中心化／缩放、群组内变化和模型设计，再将更换 optimizer 作为数值敏感性检查。warning 消失不等于科学问题已解决。

Singular fit 通常表现为随机方差接近零或相关接近 ±1，说明当前数据把某个维度估在边界。它不是机械删除全部随机斜率的命令；应结合预先问题、设计支持、不同合理结构和敏感性结果透明决定。

11 缺失、选择与因果解释的边界

11.1 Mixed model 不会自动修复缺失

最大似然 mixed model 可以使用不同时点数的个体，但无偏解释仍依赖缺失机制和模型。常见 likelihood 推断要求：给定已观测历史和模型协变量后，缺失可合理视为 MAR。

  • 症状恶化者更容易停止日记时，普通模型可能有偏；
  • 整家医院退出与单名患者失访是不同层级的选择；
  • MNAR dropout 需要 pattern-mixture、selection model 或其他敏感性分析；
  • 单次均值填补会扭曲方差与层级关系；
  • 群组大小若与潜在结局相关，可能存在 informative cluster size。

11.2 处理依赖不等于识别因果效应

多层模型可以表示相关性，却不会自动提供：

  • 清晰的干预策略与共同时间零点；
  • 条件可交换性或无未测量混杂；
  • positivity；
  • 一致性与无干扰；
  • 对时间变化治疗与时间变化混杂的正确处理；
  • 对选择进入机构、失访或测量误差的修复。

医院层干预尤其需要医院层混杂控制和足够医院重叠。患者层显著系数也不应自动翻译成患者改变暴露后的因果效应。完整框架见 因果推断详解。

随机截距不是“未测混杂控制器” 随机截距描述未观测群组偏离的分布，并不保证它与模型协变量独立，也不自动控制所有群组层共同原因。若 random effect 与协变量相关，可考虑加入群组均值的 correlated-random-effects／Mundlak 分解，并说明额外假设。

12 四个应用结果整合

12.1 动态结果摘要

case_summary <- data.frame(
  应用 = c("患者—诊所", "学生—学校", "时点—个体—治疗师", "患者—医院感染"),
  数据结构 = c(
    sprintf("%d 人／%d 诊所", nrow(med), length(unique(med$clinic))),
    sprintf("%d 人／%d 学校", nrow(soc), length(unique(soc$school))),
    sprintf(
      "%d 次／%d 人／%d 治疗师",
      nrow(psy), length(unique(psy$person)), length(unique(psy$therapist))
    ),
    sprintf("%d 人／%d 医院", nrow(bin), length(unique(bin$hospital)))
  ),
  主要结果 = c(
    sprintf("干预均值差 %.2f mmHg", med_fixed$估计值[2]),
    sprintf(
      "within SES %.2f；跨层交互 %.2f",
      soc_fixed$估计值[2], soc_fixed$估计值[5]
    ),
    sprintf("干预额外 session 斜率 %.2f", psy_fixed$估计值[6]),
    sprintf(
      "检查清单条件 OR %.2f；参考患者医院效应边际化风险差 %.3f",
      bin_results$条件OR[2], diff(p_marginal)
    )
  ),
  关键异质性 = c(
    sprintf("空模型 ICC %.3f", med_empty_icc),
    sprintf("空模型 ICC %.3f；SES slope SD %.2f", soc_empty_icc, soc_vc_table$数值[2]),
    sprintf("治疗师 SD %.2f；个体 slope SD %.2f", psy_vc_table$数值[4], psy_vc_table$数值[2]),
    sprintf("latent ICC %.3f；MOR %.2f", icc_latent, mor)
  ),
  check.names = FALSE
)

knitr::kable(case_summary, caption = "四个多层模型应用的动态结果摘要")
四个多层模型应用的动态结果摘要
应用 数据结构 主要结果 关键异质性
患者—诊所 1425 人／48 诊所 干预均值差 -4.11 mmHg 空模型 ICC 0.174
学生—学校 2098 人／70 学校 within SES 2.77；跨层交互 -1.01 空模型 ICC 0.437；SES slope SD 1.60
时点—个体—治疗师 1120 次／160 人／20 治疗师 干预额外 session 斜率 -0.40 治疗师 SD 1.64；个体 slope SD 0.29
患者—医院感染 2501 人／56 医院 检查清单条件 OR 0.57；参考患者医院效应边际化风险差 -0.048 latent ICC 0.095；MOR 1.75

这些结果使用不同 estimand 和尺度，不能按数值大小横向比较。连续模型给均值差与斜率，logistic GLMM 给条件 OR；ICC 在随机斜率和非高斯模型中也具有不同定义。

13 可审核报告模板

共分析 [N] 次观测，来自 [P] 名个体和 [J] 个群组。结局为 [定义与尺度]，时间以 [零点] 中心化。模型固定部分包含 [变量、非线性与交互]，随机部分包含 [群组截距、个体截距、随机斜率及相关]。连续结局使用 [REML／ML] 拟合，模型比较使用 [方法]；非高斯结局采用 [link 与 family]。群组、个体和残差标准差分别为 [数值]，[ICC／VPC／MOR] 为 [数值与明确定义]。主要 [均值差／斜率／条件 OR／边际风险差] 为 [估计与 95% CI]。诊断显示 [收敛、奇异性、残差、影响群组和过度离散]。结论主要推广至 [目标个体与群组范围]，并依赖 [随机效应、缺失、选择和因果假设]。

至少同时报告：

  • 每一层的单位数、群组大小范围和不平衡；
  • predictor 的层级、中心化方式与零点；
  • 完整 fixed 与 random formula；
  • random-effect 标准差、相关和残差尺度；
  • 条件还是边际效应及其尺度；
  • ML／REML、optimizer、区间和自由度方法；
  • 收敛、singularity、残差和影响群组诊断；
  • 缺失、序列相关、未测混杂和外推限制。

14 常见错误速查

common_errors <- data.frame(
  错误 = c(
    "忽略群组后把所有行当独立",
    "看到 ID 就机械添加随机截距",
    "把 fixed／random 理解为变量是否变化",
    "ICC 小就忽略聚类",
    "把 ICC 写成群组因果贡献",
    "用个体数代替群组数评价高层效应",
    "ID 只在群组内唯一却当全局唯一",
    "把 crossed 结构错写成 nested",
    "不分解时间变化变量的 within／between 效应",
    "以为总体均值中心化已经分离两个效应",
    "需要随机斜率却只放随机截距",
    "把随机斜率当成 AR(1) 残差结构",
    "比较不同 fixed effects 的 REML likelihood",
    "只凭逐步 p 值选择随机结构",
    "忽略 convergence 或 singular warning",
    "把 GLMM 条件 OR 当成风险比或边际 OR",
    "用 random effects 给机构确定排名",
    "认为 mixed model 自动处理所有缺失",
    "把随机截距当作未测混杂控制器"
  ),
  更好的做法 = c(
    "按设计表示相关性并检查每层样本量",
    "先说明共享机制和推断单位",
    "用共同系数与群组偏离分布解释",
    "结合群组大小、群组数和暴露层级",
    "称为模型方差／相关结构",
    "高层精度主要由高层单位数决定",
    "创建全局唯一 ID 或显式交互 ID",
    "分别建模学校和社区随机效应",
    "加入群组内偏离与群组均值",
    "按问题选择 group-mean decomposition",
    "有设计支持时拟合并诊断 random slope",
    "另查剩余时间序列相关",
    "用 ML 重拟合后比较固定部分",
    "由设计与科学问题预先确定候选结构",
    "排查尺度、信息和模型复杂度并透明报告",
    "报告正确尺度并计算目标边际量",
    "显示收缩、不确定性和病例结构",
    "说明 MAR 等假设并做 MNAR 敏感性分析",
    "明确 random-effect 与协变量独立等额外假设"
  ),
  check.names = FALSE
)

knitr::kable(common_errors, caption = "多层模型常见错误及更好的做法")
多层模型常见错误及更好的做法
错误 更好的做法
忽略群组后把所有行当独立 按设计表示相关性并检查每层样本量
看到 ID 就机械添加随机截距 先说明共享机制和推断单位
把 fixed／random 理解为变量是否变化 用共同系数与群组偏离分布解释
ICC 小就忽略聚类 结合群组大小、群组数和暴露层级
把 ICC 写成群组因果贡献 称为模型方差／相关结构
用个体数代替群组数评价高层效应 高层精度主要由高层单位数决定
ID 只在群组内唯一却当全局唯一 创建全局唯一 ID 或显式交互 ID
把 crossed 结构错写成 nested 分别建模学校和社区随机效应
不分解时间变化变量的 within／between 效应 加入群组内偏离与群组均值
以为总体均值中心化已经分离两个效应 按问题选择 group-mean decomposition
需要随机斜率却只放随机截距 有设计支持时拟合并诊断 random slope
把随机斜率当成 AR(1) 残差结构 另查剩余时间序列相关
比较不同 fixed effects 的 REML likelihood 用 ML 重拟合后比较固定部分
只凭逐步 p 值选择随机结构 由设计与科学问题预先确定候选结构
忽略 convergence 或 singular warning 排查尺度、信息和模型复杂度并透明报告
把 GLMM 条件 OR 当成风险比或边际 OR 报告正确尺度并计算目标边际量
用 random effects 给机构确定排名 显示收缩、不确定性和病例结构
认为 mixed model 自动处理所有缺失 说明 MAR 等假设并做 MNAR 敏感性分析
把随机截距当作未测混杂控制器 明确 random-effect 与协变量独立等额外假设

15 练习与答案

  1. 医院空模型 ICC=0.08，能否写成“医院造成 8% 的感染”？
  2. 50 家学校各 20 人与 5 家学校各 200 人，哪一个更支持学校层资源效应？
  3. 日记压力不分解就进入模型，混合了哪两个问题？
  4. 随机截距模型是否允许时间斜率因人而异？
  5. 截距–斜率相关为负，为什么必须先说明时间零点？
  6. predict(..., re.form = NA) 在 logistic GLMM 中是否就是总体边际风险？
  7. 模型出现 singular fit，是否应自动删除所有随机斜率？
  8. Mixed model 使用了所有可用随访记录，是否已经消除失访偏倚？
  9. 学生同时属于学校与居住社区且二者不嵌套，应采用什么结构？
  10. 调整后 mixed-model treatment 系数何时才可作因果解释？
显示练习答案
  1. 不能。ICC 是模型尺度的方差／相关结构，不是医院的因果贡献。
  2. 通常是 50 家学校，因为学校层独立信息更多。
  3. 同一人今天比自己通常水平更高的 within-person 关系，以及长期高压力人与低压力人的 between-person 关系。
  4. 不允许；需要个体时间随机斜率。
  5. 截距表示 time=0，改变零点会改变截距及其与斜率的相关解释。
  6. 不是。它把随机效应设为零；边际风险需要对随机效应分布积分并明确协变量分布。
  7. 不应机械删除。先检查群组内变化、尺度、样本支持、相关结构和预先科学问题。
  8. 没有。通常仍依赖给定模型信息后的 MAR 等缺失假设，MNAR 需要敏感性分析。
  9. 使用学校与社区的 crossed random effects，而非强行嵌套。
  10. 只有在干预、时间零点、可交换性、positivity、一致性、无干扰和选择机制等识别条件可信时。

16 快速参考

quick_reference <- data.frame(
  目标 = c(
    "两层随机截距 LMM", "随机斜率", "三层纵向模型",
    "二分类 GLMM", "提取固定效应", "方差成分",
    "随机效应条件模态", "奇异性", "典型群组预测", "已有群组条件预测"
  ),
  `R 入口` = c(
    "lme4::lmer(y ~ x + (1 | cluster), data=d)",
    "lme4::lmer(y ~ x + (1 + x | cluster), data=d)",
    "lme4::lmer(y ~ time*A + (1+time|person) + (1|site), data=d)",
    "lme4::glmer(y ~ x + (1|cluster), family=binomial, data=d)",
    "lme4::fixef(model)",
    "lme4::VarCorr(model)",
    "lme4::ranef(model, condVar=TRUE)",
    "lme4::isSingular(model)",
    "predict(model, re.form=NA)",
    "predict(model, re.form=NULL)"
  ),
  提醒 = c(
    "报告每层样本量和 ICC",
    "需要群组内 x 变化与足够群组",
    "ID 必须正确，检查剩余序列相关",
    "效应通常是条件 OR",
    "固定效应不是无条件边际效应的同义词",
    "报告 SD、相关和 residual scale",
    "不是无误差的群组真值",
    "边界诊断，不是自动删项命令",
    "GLMM 中不等于积分后的边际均值",
    "只适用于已有且已估计随机效应的群组"
  ),
  check.names = FALSE
)

knitr::kable(quick_reference, caption = "多层模型核心 R 入口与解释提醒")
多层模型核心 R 入口与解释提醒
目标 R 入口 提醒
两层随机截距 LMM lme4::lmer(y ~ x + (1 | cluster), data=d) 报告每层样本量和 ICC
随机斜率 lme4::lmer(y ~ x + (1 + x | cluster), data=d) 需要群组内 x 变化与足够群组
三层纵向模型 lme4::lmer(y ~ time*A + (1+time|person) + (1|site), data=d) ID 必须正确，检查剩余序列相关
二分类 GLMM lme4::glmer(y ~ x + (1|cluster), family=binomial, data=d) 效应通常是条件 OR
提取固定效应 lme4::fixef(model) 固定效应不是无条件边际效应的同义词
方差成分 lme4::VarCorr(model) 报告 SD、相关和 residual scale
随机效应条件模态 lme4::ranef(model, condVar=TRUE) 不是无误差的群组真值
奇异性 lme4::isSingular(model) 边界诊断，不是自动删项命令
典型群组预测 predict(model, re.form=NA) GLMM 中不等于积分后的边际均值
已有群组条件预测 predict(model, re.form=NULL) 只适用于已有且已估计随机效应的群组

最小分析工作流

# 1. 审计每层单位、ID、群组大小和 predictor 的层内变化。
table(d$cluster_id)
tapply(d$x, d$cluster_id, stats::var)

# 2. 先拟合空模型并分解方差。
m0 <- lme4::lmer(y ~ 1 + (1 | cluster_id), data = d, REML = TRUE)
lme4::VarCorr(m0)

# 3. 根据问题分解 within/between，并写出有依据的随机结构。
m1 <- lme4::lmer(
  y ~ x_within + x_cluster_mean + z_cluster +
    x_within:z_cluster + (1 + x_within | cluster_id),
  data = d,
  REML = TRUE,
  control = lme4::lmerControl(optimizer = "bobyqa")
)

# 4. 检查 fixed effects、variance components、singularity、残差和影响群组。
lme4::fixef(m1)
lme4::VarCorr(m1)
lme4::isSingular(m1)

# 5. 区分已有群组条件预测、典型群组预测和新群组预测。
predict(m1, re.form = NULL)
predict(m1, re.form = NA)

16.1 最终检查清单

  • 是否明确每行、每层和推断对象？
  • ID 是否全局唯一，嵌套／交叉／多重成员结构是否正确？
  • 是否报告每层样本量、群组大小范围和不平衡？
  • predictor 是否按科学问题分成 within 与 between 成分？
  • 零点和中心化方式是否使截距与交互可解释？
  • random slope 是否有群组内变化与群组数支持？
  • ML／REML、optimizer、区间和自由度方法是否预先说明？
  • 是否检查 convergence、singularity、残差、序列相关和影响群组？
  • GLMM 是否区分条件、典型群组与边际结果？
  • 随机效应是否连同收缩和不确定性解释，而非直接排名？
  • 缺失、信息性群组大小、测量误差与选择是否另行评估？
  • 因果语言是否与研究设计和识别假设一致？

下一步学习方向

后续可学习非线性时间轨迹、样条 random slopes、异方差与 AR(1) 残差、cross-classified 与 multiple-membership 模型、Bayesian multilevel models、joint longitudinal–survival models、multilevel mediation、cluster randomized trials、小样本自由度校正、parametric bootstrap、survey-weighted multilevel models，以及模型的外部验证与动态预测。

sessionInfo()
## 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] nlme_3.1-169     cli_3.6.6        knitr_1.51       rlang_1.3.0
##  [5] xfun_0.60        reformulas_0.4.4 jsonlite_2.0.0   minqa_1.2.8
##  [9] htmltools_0.5.9  lme4_2.0-6       sass_0.4.10      rmarkdown_2.31
## [13] grid_4.6.1       evaluate_1.0.5   jquerylib_0.1.4  MASS_7.3-65
## [17] fastmap_1.2.0    yaml_2.3.12      lifecycle_1.0.5  compiler_4.6.1
## [21] Rcpp_1.1.2       lattice_0.22-9   digest_0.6.39    nloptr_2.2.1
## [25] R6_2.6.1         Rdpack_2.6.6     splines_4.6.1    rbibutils_2.4.1
## [29] bslib_0.12.0     Matrix_1.7-5     tools_4.6.1      boot_1.3-32
## [33] cachem_1.1.0

教学用途：本模块用于学习多层数据的建模、诊断、解释与报告，不能替代具体研究的统计分析计划、数据审计或专业统计咨询。所有数据和数值结果均为模拟结果。