适用对象公共卫生、医学、实验室科学与健康服务研究学习者

学习时长约 180–240 分钟

先修要求描述统计、置信区间、线性模型与基本 R 语法

贯穿案例一家公共卫生机构在 12 家诊所检验两种可组合的降压支持：健康教练（coaching）和居家血压监测（home monitoring）。每家诊所将 16 名参与者随机分到四个组合，每组合 4 人，共 192 人。响应是 12 周收缩压下降值（mmHg，越大越好）。数据由固定种子模拟，不是临床资料；在生成模型的 0/1 编码下，无另一干预时教练和监测的简单效应分别为 3.0 和 2.2 mmHg，差异之差交互作用为 2.8 mmHg，另有诊所差异和个体误差。因此生成模型中的平均主效应分别为 4.4 和 3.6 mmHg，须由四个 cell mean 计算。

随机化不是一句标签只有分配机制真正随机、分配隐藏得到执行、分析保留随机分组且干预间污染受到控制时，随机实验的因果解释才有基础。事后把观察性资料称为“自然实验”，或仅在回归中加入协变量，都不能复现随机分配。

如何使用本教程

建议按下面的研究链条学习：

问题与 estimand → 实验单位 → 随机化/重复/局部控制 → CRD 与 RCBD → 2×2 析因设计 → 估计和比较 → 诊断与随机化推断 → 功效 → 实施偏差 → 更复杂设计

学习目标

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

区分实验单位、观察单位、处理、因子、水平、响应和 estimand；
用随机化、重复和局部控制解释实验为何能够支持有效比较；
为完全随机设计（CRD）和随机完全区组设计（RCBD）生成可审核的分配表；
解释 $2\times2$ 析因设计中的两个主效应和差异之差交互作用；
用 lm()、平衡设计 ANOVA、模型矩阵和线性组合估计四个 cell mean 及置信区间；
区分预设对比与 Tukey 事后多重比较，避免按 P 值挑结果；
比较区组设计与忽略区组分析的精度，并检查模型残差和影响点；
实施保持诊所内分配数目的 restricted randomization permutation test；
用解析近似和模拟评估样本量/功效，同时透明陈述假设；
按 ITT 原则处理不依从，识别缺失、污染和伪重复造成的偏差；
判断何时需要 cluster、split-plot、crossover、Latin square 或重复测量设计；
为连续、二元和计数结局选择与分配单位一致的分析方法。

1 从研究问题到可估计效应

1.1 DOE 的核心术语

实验设计（design of experiments, DOE）是在观察响应之前安排处理、随机化、重复和测量的规则。它的目的不是让数据“看起来平衡”，而是让比较具有已知的分配机制和可解释的不确定性。

术语	本案例中的定义	审核问题
实验单位	被独立随机分配的一名参与者	干预实际上能否独立施加到每个人？
观察单位	12 周时的一次个体结局记录	是否误把一人的重复记录当成独立实验单位？
处理	教练与居家监测的一个组合	对照组获得什么常规照护？
因子与水平	两个因子，各为 No/Yes 两水平	水平是否可实施、可区分且保持一致？
区组	诊所	区组在处理前形成，且区组内更同质吗？
响应	12 周收缩压下降 mmHg	时间点、测量姿势和设备是否预先规定？
estimand	某个目标总体、依从策略和比较下的平均效果	比较谁、在何种条件、用什么尺度？

本页的主要 estimand 是：在参与研究的 12 家诊所等权平均，按最初随机分配（ITT），教练和监测对 12 周血压下降均值的主效应及其交互作用。若目标是推广到全国所有诊所，就需要把诊所视为来自目标总体的样本，并改变设计或模型；12 个固定诊所本身不会自动支持该推广。

1.2 四个因果对比

记 $\mu_{cm}=E\{Y(c,m)\}$ ，其中 $c,m\in\{0,1\}$ 分别表示教练和监测。

$\begin{aligned} \text{教练平均主效应} &=\tfrac12[(\mu_{10}-\mu_{00})+(\mu_{11}-\mu_{01})],\\ \text{监测平均主效应} &=\tfrac12[(\mu_{01}-\mu_{00})+(\mu_{11}-\mu_{10})],\\ \text{交互作用} &=(\mu_{11}-\mu_{01})-(\mu_{10}-\mu_{00}). \end{aligned}$

交互作用是“一个因素的效果是否随另一个因素改变”。它不是两条线肉眼是否相交的标签，也不是“联合组显著而单独组不显著”。本案例设定正交互 2.8 mmHg：监测开启时，教练效果比监测关闭时额外增加 2.8 mmHg。

2 六个设计支柱

2.1 随机化、重复与局部控制

随机化：让未控制预后因素在分配机制下可比较，并为随机化推断提供参考分布。
独立重复：每个处理组合需要多个真正独立的实验单位，以估计个体变异。技术重复提高测量精度，但不增加实验单位数。
局部控制：先按诊所等强预后变量形成区组，再在区组内随机化，减少残差变异。
对照：明确常规照护、安慰处理或活性对照，避免将额外关注与干预成分混为一谈。
盲法与测量标准化：参与者未必能对行为干预盲法，但结局测量者、实验室人员和分析标签可尽量盲法。
预注册与分配隐藏：在看结局前固定主要 estimand、结局、分析集、对比和缺失策略；使用不可预测的中央分配或封闭流程。

随机化清单要保存什么？保存随机种子只是最低要求。还应保存代码版本、区组形成时间、分配比例、限制条件、随机列表生成者、分配执行者、揭盲时间、偏离原因与审计日志。生产环境中不要把公开教学种子当成不可预测的分配工具。

2.2 完全随机设计（CRD）

CRD 把所有实验单位放入一个随机化池，适用于实验单位较同质、没有强预后分层因素且处理可以独立施加的场景。若 20 人平均分到四组，可以这样生成一个可复现的教学列表：

set.seed(20260913)
crd_example <- data.frame(
  入组顺序 = sprintf("S%02d", seq_len(20)),
  处理组合 = sample(rep(c("C0_M0", "C1_M0", "C0_M1", "C1_M1"), each = 5))
)
knitr::kable(crd_example, caption = "20 人四组等比例 CRD 教学分配表")

20 人四组等比例 CRD 教学分配表
入组顺序	处理组合
S01	C1_M1
S02	C0_M0
S03	C1_M0
S04	C0_M0
S05	C0_M1
S06	C1_M1
S07	C1_M1
S08	C1_M0
S09	C0_M0
S10	C1_M0
S11	C0_M0
S12	C0_M1
S13	C1_M1
S14	C1_M1
S15	C1_M0
S16	C0_M0
S17	C0_M1
S18	C1_M0
S19	C0_M1
S20	C0_M1

等样本数可提高比较精度，但“恰好平衡”不是随机化的定义。若流失、排除或实施失败在分配后发生，不能通过删除参与者来恢复外观平衡；应保留 ITT 分组并说明偏离。

2.3 为什么本案例采用 RCBD

诊所之间可能在服务流程、基线风险和测量团队上不同。RCBD 先把每家诊所作为一个区组，再在每个诊所内对四个组合各分 4 人。这样每个处理都在每家诊所出现，处理比较不会与诊所构成混杂。

randomization_preview <- doe_data[seq_len(32), c(
  "clinic", "participant_in_clinic", "id", "treatment"
)]
names(randomization_preview) <- c("诊所", "诊所内入组顺序", "研究编号", "处理编码")
knitr::kable(
  randomization_preview,
  caption = "前两家诊所的随机分配表：顺序随机，但每个组合各 4 人"
)

前两家诊所的随机分配表：顺序随机，但每个组合各 4 人
诊所	诊所内入组顺序	研究编号	处理编码
诊所01	1	P001	C1_M1
诊所01	2	P002	C0_M1
诊所01	3	P003	C1_M0
诊所01	4	P004	C0_M1
诊所01	5	P005	C0_M0
诊所01	6	P006	C0_M1
诊所01	7	P007	C1_M0
诊所01	8	P008	C1_M1
诊所01	9	P009	C1_M0
诊所01	10	P010	C1_M0
诊所01	11	P011	C0_M0
诊所01	12	P012	C0_M0
诊所01	13	P013	C1_M1
诊所01	14	P014	C1_M1
诊所01	15	P015	C0_M0
诊所01	16	P016	C0_M1
诊所02	1	P017	C1_M0
诊所02	2	P018	C1_M0
诊所02	3	P019	C1_M0
诊所02	4	P020	C1_M1
诊所02	5	P021	C0_M1
诊所02	6	P022	C1_M0
诊所02	7	P023	C1_M1
诊所02	8	P024	C0_M0
诊所02	9	P025	C1_M1
诊所02	10	P026	C0_M1
诊所02	11	P027	C0_M0
诊所02	12	P028	C0_M0
诊所02	13	P029	C0_M0
诊所02	14	P030	C1_M1
诊所02	15	P031	C0_M1
诊所02	16	P032	C0_M1

balance_table <- as.data.frame.matrix(xtabs(~ clinic + treatment, data = doe_data))
balance_table <- cbind(诊所 = rownames(balance_table), balance_table)
rownames(balance_table) <- NULL
knitr::kable(balance_table, caption = "每家诊所内四个处理组合的样本数审计")

每家诊所内四个处理组合的样本数审计
诊所	C0_M0	C1_M0	C0_M1	C1_M1
诊所01	4	4	4	4
诊所02	4	4	4	4
诊所03	4	4	4	4
诊所04	4	4	4	4
诊所05	4	4	4	4
诊所06	4	4	4	4
诊所07	4	4	4	4
诊所08	4	4	4	4
诊所09	4	4	4	4
诊所10	4	4	4	4
诊所11	4	4	4	4
诊所12	4	4	4	4

这种设计是受限随机化：允许的分配必须满足每诊所每组合 4 人。分析和随机化检验都要尊重同一限制。若跨全体 192 人任意置换标签，会产生实验中从未可能出现的分配。

3 2×2 析因设计：一次实验回答三个问题

3.1 为什么不做两个独立实验

四组析因设计同时估计教练主效应、监测主效应和交互作用，并让每个因素的主效应借用另一个因素两个水平的数据。若科学上合理地假设处理能共同实施，它通常比两个互不相干的试验更有效率。

但“主效应”是对另一个因素水平的平均。若交互很大，单一平均主效应会掩盖有意义的条件效果；这时应同时报告四个 cell mean 和简单效应。

3.2 数据字典与生成机制

data_dictionary <- data.frame(
  变量 = c("clinic", "id", "coaching", "home_monitoring", "treatment", "bp_reduction"),
  含义 = c(
    "随机化区组：12 家诊所", "个体研究编号", "是否分配健康教练",
    "是否分配居家监测", "四个析因组合", "12 周收缩压下降 mmHg；越大越好"
  ),
  角色 = c("设计变量", "实验单位标识", "因素 A", "因素 B", "处理", "连续响应")
)
knitr::kable(data_dictionary, caption = "共享主案例的数据字典")

共享主案例的数据字典
变量	含义	角色
clinic	随机化区组：12 家诊所	设计变量
id	个体研究编号	实验单位标识
coaching	是否分配健康教练	因素 A
home_monitoring	是否分配居家监测	因素 B
treatment	四个析因组合	处理
bp_reduction	12 周收缩压下降 mmHg；越大越好	连续响应

模拟模型为

$Y_{ijk}=4+3C_{ijk}+2.2M_{ijk}+2.8C_{ijk}M_{ijk}+b_j+\varepsilon_{ijk},$ $b_j\sim N(0,2.2^2),\qquad \varepsilon_{ijk}\sim N(0,4.5^2).$

这里 $b_j$ 只是生成 12 家诊所差异；主分析把这些已入组诊所作为固定区组。0/1 编码下的 3.0 和 2.2 是另一干预关闭时的简单效应；对另一因素两个水平等权平均后，真实平均主效应为 $3.0+2.8/2=4.4$ 和 $2.2+2.8/2=3.6$ mmHg。响应中的随机噪声会使样本估计不等于生成参数，这正是不确定性分析存在的原因。

4 探索性分析：先看 cell mean，再看模型

4.1 四组样本量、均值与标准误

cell_split <- split(doe_data[["bp_reduction"]], doe_data[["treatment"]])
cell_summary <- data.frame(
  处理 = names(cell_split),
  n = vapply(cell_split, length, integer(1)),
  平均下降 = vapply(cell_split, mean, numeric(1)),
  标准差 = vapply(cell_split, sd, numeric(1)),
  标准误 = vapply(cell_split, function(x) sd(x) / sqrt(length(x)), numeric(1))
)
rownames(cell_summary) <- NULL
cell_summary[["下限"]] <- cell_summary[["平均下降"]] -
  qt(0.975, cell_summary[["n"]] - 1) * cell_summary[["标准误"]]
cell_summary[["上限"]] <- cell_summary[["平均下降"]] +
  qt(0.975, cell_summary[["n"]] - 1) * cell_summary[["标准误"]]
knitr::kable(cell_summary, digits = 2, caption = "四个处理组合的未调整描述统计")

四个处理组合的未调整描述统计
处理	n	平均下降	标准差	标准误	下限	上限
C0_M0	48	2.96	4.72	0.68	1.59	4.33
C1_M0	48	5.77	6.23	0.90	3.96	7.58
C0_M1	48	6.18	5.55	0.80	4.56	7.79
C1_M1	48	11.24	5.08	0.73	9.76	12.71

set_doe_plot_font()
plot(
  seq_len(nrow(cell_summary)), cell_summary[["平均下降"]],
  ylim = range(cell_summary[, c("下限", "上限")]), pch = 19,
  col = palette_doe["navy"], xaxt = "n", xlab = "处理组合",
  ylab = "12 周收缩压下降（mmHg）",
  main = "先展示绝对均值，再解释对比"
)
axis(1, at = seq_len(nrow(cell_summary)), labels = cell_summary[["处理"]])
arrows(
  seq_len(nrow(cell_summary)), cell_summary[["下限"]],
  seq_len(nrow(cell_summary)), cell_summary[["上限"]],
  angle = 90, code = 3, length = 0.06, col = palette_doe["teal"]
)
abline(h = 0, lty = 3, col = palette_doe["gray"])

点和误差线图比较无教练无监测、有教练无监测、无教练有监测和两者均有四组的平均收缩压下降值；误差线是忽略诊所区组结构的朴素描述区间，数值越高表示下降越多。

四个处理组合的未调整平均收缩压下降值与朴素描述区间

这些是未调整均值，误差线由每组总体标准差除以 $\sqrt{48}$ 得到，只是忽略诊所区组结构的朴素描述区间，不是主模型的正式置信区间。因为设计在每家诊所完全平衡，调整诊所后的处理对比会与未调整对比非常接近，但正式区组模型的标准误通常更小；后文给出与 estimand 对齐的模型矩阵区间。误差线重叠与否不是交互作用检验。

4.2 交互作用的几何解释

set_doe_plot_font()
interaction.plot(
  x.factor = doe_data[["coaching"]],
  trace.factor = doe_data[["home_monitoring"]],
  response = doe_data[["bp_reduction"]],
  fun = mean, type = "b", pch = c(16, 17), lwd = 2,
  col = c(palette_doe["orange"], palette_doe["teal"]),
  ylim = range(cell_summary[["平均下降"]]) + c(-1.2, 2.5),
  xlab = "健康教练", ylab = "平均收缩压下降（mmHg）",
  trace.label = "居家监测", fixed = TRUE, legend = FALSE,
  main = "交互作用是斜率之差，而非显著性标签之差"
)
legend(
  "topleft", legend = c("无监测", "有监测"),
  col = c(palette_doe["orange"], palette_doe["teal"]),
  pch = c(16, 17), lty = 1, lwd = 2, bty = "n", title = "居家监测"
)

交互作用折线图显示在无居家监测和有居家监测条件下，健康教练从 No 到 Yes 的平均收缩压下降变化；两条不平行的线提示效果修饰。

健康教练与居家监测对平均收缩压下降的交互作用图

若两条线平行，样本中的交互对比接近 0；若明显不平行，表示教练的简单效应随监测状态变化。图是估计与沟通工具，正式不确定性仍应来自预设的交互对比。

5 RCBD 线性模型与平衡设计 ANOVA

5.1 拟合预先规定的模型

rcbd_fit <- lm(
  bp_reduction ~ clinic + coaching * home_monitoring,
  data = doe_data
)
summary(rcbd_fit)

## 
## Call:
## lm(formula = bp_reduction ~ clinic + coaching * home_monitoring, 
##     data = doe_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.767  -3.148   0.204   3.228  12.736 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       0.320      1.315    0.24   0.8081    
## clinic诊所02                      7.536      1.663    4.53  1.1e-05 ***
## clinic诊所03                     -2.885      1.663   -1.73   0.0845 .  
## clinic诊所04                      2.694      1.663    1.62   0.1071    
## clinic诊所05                      1.904      1.663    1.14   0.2538    
## clinic诊所06                      6.850      1.663    4.12  5.8e-05 ***
## clinic诊所07                      0.191      1.663    0.12   0.9085    
## clinic诊所08                      4.397      1.663    2.64   0.0089 ** 
## clinic诊所09                      4.189      1.663    2.52   0.0127 *  
## clinic诊所10                     -0.229      1.663   -0.14   0.8906    
## clinic诊所11                      3.667      1.663    2.21   0.0287 *  
## clinic诊所12                      3.400      1.663    2.04   0.0424 *  
## coachingYes                       2.806      0.960    2.92   0.0039 ** 
## home_monitoringYes                3.214      0.960    3.35   0.0010 ***
## coachingYes:home_monitoringYes    2.255      1.358    1.66   0.0986 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.7 on 177 degrees of freedom
## Multiple R-squared:  0.459,  Adjusted R-squared:  0.416 
## F-statistic: 10.7 on 14 and 177 DF,  p-value: <2e-16

模型中的 clinic 控制区组均值；coaching * home_monitoring 展开为两个主项和一个乘积项。在默认 treatment coding 下：

coachingYes 是无监测时教练的简单效应；
home_monitoringYes 是无教练时监测的简单效应；
coachingYes:home_monitoringYes 是差异之差；
平均主效应需用明确线性组合计算，不能在有交互时直接把简单效应系数改名为“总体主效应”。

5.2 平衡设计的 ANOVA 分解

anova_table <- anova(rcbd_fit)
anova_display <- data.frame(
  来源 = rownames(anova_table),
  自由度 = anova_table[, "Df"],
  平方和 = anova_table[, "Sum Sq"],
  均方 = anova_table[, "Mean Sq"],
  F值 = ifelse(
    is.na(anova_table[, "F value"]), "—",
    formatC(anova_table[, "F value"], format = "f", digits = 3)
  ),
  P值 = ifelse(
    is.na(anova_table[, "Pr(>F)"]), "—",
    format_p(anova_table[, "Pr(>F)"])
  )
)
rownames(anova_display) <- NULL
knitr::kable(
  anova_display, digits = 3,
  align = c("l", "r", "r", "r", "r", "l"),
  caption = "RCBD 平衡析因模型的顺序 ANOVA 表"
)

RCBD 平衡析因模型的顺序 ANOVA 表
来源	自由度	平方和	均方	F值	P值
clinic	11	1617	147.0	6.645	<0.001
coaching	1	743	742.7	33.567	<0.001
home_monitoring	1	904	904.5	40.875	<0.001
coaching:home_monitoring	1	61	61.0	2.757	0.0986
Residuals	177	3917	22.1	—	—

由于每个诊所内四个组合等频，处理项与诊所正交，两个因素和交互也正交；在这个特定平衡模型中，顺序平方和不依赖这些处理项的先后顺序。失访或不等比例会破坏这种性质；此时应回到 estimand、模型矩阵和预设对比，不要机械套用“Type I/II/III”标签。

P 值不是效应是否存在的开关较大的 P 值表示数据在特定模型下没有提供足够精确的反对零假设证据；它不等于证明无效，也不证明两个效果相等。应报告估计、置信区间、临床相关尺度和设计功效。等效性或非劣效性需要预设界值与专门设计。

6 用模型矩阵得到四个调整后 cell mean

6.1 为什么手动线性组合值得学习

模型系数随参考水平和编码改变，但四个 cell mean 及其科学对比不应改变。下面对每个处理组合，在 12 家诊所分别构造模型行，再对诊所等权平均；这与本页“12 家参与诊所平均效果”的 estimand 对齐。

cell_key <- data.frame(
  coaching = factor(c("No", "Yes", "No", "Yes"), levels = levels(doe_data[["coaching"]])),
  home_monitoring = factor(c("No", "No", "Yes", "Yes"), levels = levels(doe_data[["home_monitoring"]])),
  处理 = c("C0_M0", "C1_M0", "C0_M1", "C1_M1")
)

prediction_grid <- do.call(rbind, lapply(seq_len(nrow(cell_key)), function(k) {
  data.frame(
    clinic = factor(clinic_labels, levels = levels(doe_data[["clinic"]])),
    coaching = factor(rep(as.character(cell_key[["coaching"]][k]), n_clinics),
                      levels = levels(doe_data[["coaching"]])),
    home_monitoring = factor(rep(as.character(cell_key[["home_monitoring"]][k]), n_clinics),
                             levels = levels(doe_data[["home_monitoring"]])),
    cell_id = k
  )
}))

X_grid <- model.matrix(
  delete.response(terms(rcbd_fit)), prediction_grid,
  contrasts.arg = rcbd_fit[["contrasts"]]
)
L_cells <- t(vapply(seq_len(nrow(cell_key)), function(k) {
  colMeans(X_grid[prediction_grid[["cell_id"]] == k, , drop = FALSE])
}, numeric(ncol(X_grid))))
colnames(L_cells) <- colnames(X_grid)

cell_estimate <- as.vector(L_cells %*% coef(rcbd_fit))
cell_se <- sqrt(diag(L_cells %*% vcov(rcbd_fit) %*% t(L_cells)))
cell_critical <- qt(0.975, df.residual(rcbd_fit))
adjusted_cells <- data.frame(
  处理 = cell_key[["处理"]],
  调整后均值 = cell_estimate,
  标准误 = cell_se,
  下限 = cell_estimate - cell_critical * cell_se,
  上限 = cell_estimate + cell_critical * cell_se
)
knitr::kable(
  adjusted_cells, digits = 2,
  caption = "对 12 家诊所等权平均的四个模型调整后 cell mean 与 95% 置信区间"
)

对 12 家诊所等权平均的四个模型调整后 cell mean 与 95% 置信区间
处理	调整后均值	标准误	下限	上限
C0_M0	2.96	0.68	1.62	4.30
C1_M0	5.77	0.68	4.43	7.11
C0_M1	6.18	0.68	4.84	7.52
C1_M1	11.24	0.68	9.90	12.58

这些区间是各 cell mean 的逐项区间，不是所有四个均值的同时置信带。若报告重点是预设对比，应直接为对比构造区间；若要探索全部两两比较，则应控制多重性。

6.2 预设主效应与交互对比

contrast_weights <- rbind(
  `教练平均主效应` = c(-0.5, 0.5, -0.5, 0.5),
  `监测平均主效应` = c(-0.5, -0.5, 0.5, 0.5),
  `交互作用（差异之差）` = c(1, -1, -1, 1),
  `联合干预 vs 均无` = c(-1, 0, 0, 1)
)
L_contrasts <- contrast_weights %*% L_cells
contrast_estimate <- as.vector(L_contrasts %*% coef(rcbd_fit))
contrast_se <- sqrt(diag(L_contrasts %*% vcov(rcbd_fit) %*% t(L_contrasts)))
contrast_t <- contrast_estimate / contrast_se
planned_results <- data.frame(
  对比 = rownames(contrast_weights),
  估计_mmHg = contrast_estimate,
  标准误 = contrast_se,
  下限 = contrast_estimate - cell_critical * contrast_se,
  上限 = contrast_estimate + cell_critical * contrast_se,
  名义未调整P值 = format_p(
    2 * pt(abs(contrast_t), df.residual(rcbd_fit), lower.tail = FALSE)
  )
)
rownames(planned_results) <- NULL
knitr::kable(
  planned_results, digits = 3,
  caption = "预先规定的析因线性对比；P 值为逐项名义值，须结合预设多重性策略"
)

预先规定的析因线性对比；P 值为逐项名义值，须结合预设多重性策略
对比	估计_mmHg	标准误	下限	上限	名义未调整P值
教练平均主效应	3.93	0.679	2.594	5.27	<0.001
监测平均主效应	4.34	0.679	3.001	5.68	<0.001
交互作用（差异之差）	2.25	1.358	-0.425	4.93	0.0986
联合干预 vs 均无	8.28	0.960	6.380	10.17	<0.001

预设对比应由科学问题决定。若“联合干预 vs 均无”是唯一主要比较，可为它分配主要一类错误率；三个析因效应可作为共同主要假设或有层级的次要假设。无论策略如何，都应在看结果前写明。

6.3 Tukey 全部两两比较何时适用

若目标确实是探索四个处理组合的所有六个两两差异，可在包含区组的 aov 模型上用 Tukey 家族区间：

tukey_fit <- aov(bp_reduction ~ clinic + treatment, data = doe_data)
tukey_matrix <- TukeyHSD(tukey_fit, which = "treatment")[["treatment"]]
tukey_results <- data.frame(
  比较 = rownames(tukey_matrix),
  差值 = tukey_matrix[, "diff"],
  下限 = tukey_matrix[, "lwr"],
  上限 = tukey_matrix[, "upr"],
  调整后P值 = format_p(tukey_matrix[, "p adj"])
)
rownames(tukey_results) <- NULL
knitr::kable(tukey_results, digits = 3, caption = "四个处理组合的 Tukey 同时比较")

四个处理组合的 Tukey 同时比较
比较	差值	下限	上限	调整后P值
C1_M0-C0_M0	2.806	0.316	5.30	0.02031
C0_M1-C0_M0	3.214	0.723	5.70	0.00547
C1_M1-C0_M0	8.275	5.784	10.77	< 0.001
C0_M1-C1_M0	0.407	-2.083	2.90	0.97430
C1_M1-C1_M0	5.468	2.978	7.96	< 0.001
C1_M1-C0_M1	5.061	2.571	7.55	< 0.001

Tukey 回答所有 cell mean 两两比较，不替代析因主效应和交互的科学解释。不要先看六个未调整 P 值，再只报告最小者；那会使名义错误率失真。

7 区组如何提高精度

7.1 与忽略诊所的 CRD 分析比较

因为处理在每家诊所平衡，忽略诊所通常不会改变处理点估计，但会把可解释的诊所差异留在残差中。下面比较教练平均主效应的标准误和残差均方。

crd_fit <- lm(bp_reduction ~ coaching * home_monitoring, data = doe_data)
crd_grid <- cell_key[, c("coaching", "home_monitoring")]
L_cells_crd <- model.matrix(delete.response(terms(crd_fit)), crd_grid,
                            contrasts.arg = crd_fit[["contrasts"]])

coach_weights <- contrast_weights["教练平均主效应", , drop = FALSE]
L_coach_rcbd <- coach_weights %*% L_cells
L_coach_crd <- coach_weights %*% L_cells_crd
precision_comparison <- data.frame(
  分析 = c("忽略诊所（当作 CRD）", "控制诊所（RCBD）"),
  教练平均效应 = c(
    as.vector(L_coach_crd %*% coef(crd_fit)),
    as.vector(L_coach_rcbd %*% coef(rcbd_fit))
  ),
  标准误 = c(
    sqrt(L_coach_crd %*% vcov(crd_fit) %*% t(L_coach_crd)),
    sqrt(L_coach_rcbd %*% vcov(rcbd_fit) %*% t(L_coach_rcbd))
  ),
  残差均方 = c(
    deviance(crd_fit) / df.residual(crd_fit),
    deviance(rcbd_fit) / df.residual(rcbd_fit)
  ),
  残差自由度 = c(df.residual(crd_fit), df.residual(rcbd_fit))
)
knitr::kable(precision_comparison, digits = 3, caption = "控制区组前后的精度比较")

控制区组前后的精度比较
分析	教练平均效应	标准误	残差均方	残差自由度
忽略诊所（当作 CRD）	3.93	0.783	29.4	188
控制诊所（RCBD）	3.93	0.679	22.1	177

区组也会消耗自由度；若区组与响应几乎无关，收益可能很小。区组变量应在处理前定义，通常选强预后因素，并避免把每个极小组合都切成区组而导致实施困难。

8 模型诊断：设计正确仍要检查测量与模型

8.1 残差、正态性、方差与影响

线性模型对均值结构、独立误差、近似恒定方差和用于小样本推断的正态误差有要求。随机化保护处理比较免受系统性基线混杂，但不会自动修复错误记录、重尾误差、地板/天花板效应或处理依赖的方差。

set_doe_plot_font()
old_par <- par(mfrow = c(2, 2), mar = c(4, 4, 2.3, 1))
plot(fitted(rcbd_fit), residuals(rcbd_fit), pch = 16,
     col = grDevices::adjustcolor(palette_doe["navy"], 0.55),
     xlab = "拟合值", ylab = "残差", main = "残差 vs 拟合值")
abline(h = 0, lty = 2, col = palette_doe["vermillion"])
qqnorm(rstandard(rcbd_fit), pch = 16,
       col = grDevices::adjustcolor(palette_doe["teal"], 0.6),
       main = "标准化残差 Q-Q")
qqline(rstandard(rcbd_fit), col = palette_doe["vermillion"], lwd = 2)
plot(fitted(rcbd_fit), sqrt(abs(rstandard(rcbd_fit))), pch = 16,
     col = grDevices::adjustcolor(palette_doe["orange"], 0.6),
     xlab = "拟合值", ylab = expression(sqrt("|标准化残差|")),
     main = "尺度—位置")
cook <- cooks.distance(rcbd_fit)
plot(cook, type = "h", col = palette_doe["purple"],
     xlab = "观察编号", ylab = "Cook 距离", main = "影响诊断")
abline(h = 4 / nrow(doe_data), lty = 2, col = palette_doe["vermillion"])

四联图依次展示拟合值对残差、残差正态 Q-Q 图、拟合值对平方根标准化残差以及各观察的 Cook 距离，用于检查均值结构、尾部、异方差和高影响观察。

RCBD 析因线性模型的残差与影响诊断

par(old_par)

4/n 只是筛查线，不是自动删除规则。发现异常点时应回到原始记录、测量设备和方案偏离；同时报告包含与不包含合理修正的敏感性分析。不能因为结果变得“不显著”或“显著”而决定删点。

diagnostic_audit <- data.frame(
  指标 = c("残差标准差", "最大绝对标准化残差", "超过 4/n 的 Cook 距离数", "最大杠杆值"),
  数值 = c(
    sigma(rcbd_fit), max(abs(rstandard(rcbd_fit))),
    sum(cooks.distance(rcbd_fit) > 4 / nrow(doe_data)),
    max(hatvalues(rcbd_fit))
  )
)
knitr::kable(diagnostic_audit, digits = 3, caption = "模型诊断的数值审计摘要")

模型诊断的数值审计摘要
指标	数值
残差标准差	4.704
最大绝对标准化残差	2.820
超过 4/n 的 Cook 距离数	10.000
最大杠杆值	0.078

9 Restricted randomization permutation test

9.1 用真实允许的分配构造零分布

随机化检验可直接利用分配机制。这里以交互作用差异之差为预设统计量，检验“对每名参与者，四个组合的潜在结局完全相同”的 Fisher sharp global null。它不是“总体平均交互作用为 0”的弱零假设检验。每次只在同一家诊所内重排标签，因此保持每组合 4 人。

interaction_statistic <- function(labels, outcome) {
  means <- tapply(outcome, factor(labels, levels = 1:4), mean)
  unname((means[4] - means[3]) - (means[2] - means[1]))
}

omnibus_treatment_f <- function(labels, outcome, clinics) {
  test_data <- data.frame(
    outcome = outcome,
    clinic = clinics,
    treatment = factor(labels, levels = 1:4)
  )
  reduced_fit <- lm(outcome ~ clinic, data = test_data)
  full_fit <- lm(outcome ~ clinic + treatment, data = test_data)
  unname(anova(reduced_fit, full_fit)[["F"]][2])
}

observed_interaction <- interaction_statistic(doe_data[["combination"]], doe_data[["bp_reduction"]])
observed_omnibus_f <- omnibus_treatment_f(
  doe_data[["combination"]], doe_data[["bp_reduction"]], doe_data[["clinic"]]
)
clinic_indices <- split(seq_len(nrow(doe_data)), doe_data[["clinic"]])

set.seed(20260915)
B <- 999
permuted_interaction <- numeric(B)
permuted_omnibus_f <- numeric(B)
for (b in seq_len(B)) {
  permuted_labels <- integer(nrow(doe_data))
  for (idx in clinic_indices) {
    permuted_labels[idx] <- sample(doe_data[["combination"]][idx])
  }
  permuted_interaction[b] <- interaction_statistic(permuted_labels, doe_data[["bp_reduction"]])
  permuted_omnibus_f[b] <- omnibus_treatment_f(
    permuted_labels, doe_data[["bp_reduction"]], doe_data[["clinic"]]
  )
}
interaction_randomization_p <-
  (1 + sum(abs(permuted_interaction) >= abs(observed_interaction))) / (B + 1)
omnibus_randomization_p <-
  (1 + sum(permuted_omnibus_f >= observed_omnibus_f)) / (B + 1)

permutation_result <- data.frame(
  预设统计量 = c("交互作用差异之差（双侧）", "诊所调整的 3-df 总体处理 F"),
  观察值 = c(observed_interaction, observed_omnibus_f),
  重随机化次数 = c(B, B),
  Fisher_sharp_null_P值 = c(interaction_randomization_p, omnibus_randomization_p)
)
knitr::kable(
  permutation_result, digits = 3,
  caption = "同一 Fisher sharp global null 下两种预设统计量的诊所内受限随机化检验"
)

同一 Fisher sharp global null 下两种预设统计量的诊所内受限随机化检验
预设统计量	观察值	重随机化次数	Fisher_sharp_null_P值
交互作用差异之差（双侧）	2.25	999	0.158
诊所调整的 3-df 总体处理 F	25.73	999	0.001

set_doe_plot_font()
old_par <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(
  permuted_interaction, breaks = 28,
  col = grDevices::adjustcolor(palette_doe["sky"], 0.75),
  border = "white", xlab = "置换交互作用（mmHg）",
  main = "零分布必须复刻实际随机化限制"
)
abline(v = c(-abs(observed_interaction), abs(observed_interaction)),
       col = palette_doe["vermillion"], lwd = 2, lty = 2)
hist(
  permuted_omnibus_f, breaks = 28,
  col = grDevices::adjustcolor(palette_doe["green"], 0.65),
  border = "white", xlab = "置换总体 F",
  main = "3-df 总体处理 F"
)
usr <- par("usr")
arrows(
  x0 = usr[2] - 0.04 * diff(usr[1:2]), y0 = usr[3] + 0.12 * diff(usr[3:4]),
  x1 = usr[2], y1 = usr[3] + 0.12 * diff(usr[3:4]),
  col = palette_doe["vermillion"], lwd = 2, length = 0.08
)
text(
  usr[2] - 0.04 * diff(usr[1:2]), usr[3] + 0.22 * diff(usr[3:4]),
  labels = sprintf("观察 F = %.2f（图外）", observed_omnibus_f),
  col = palette_doe["vermillion"], adj = 1, cex = 0.82
)

并排两个直方图分别显示九百九十九次诊所内受限重随机化的交互作用差异之差和诊所调整三自由度总体处理 F；红线标记各自观察统计量。

同一诊所内受限重随机化产生的交互作用与总体处理 F 零分布

par(old_par)

加 1 的 P 值公式避免在有限 Monte Carlo 抽样中报告 0。两行随机化检验使用同一个 Fisher sharp global null 和同一批允许分配，但统计量不同：差异之差聚焦加性交互，总体 F 对任何四组均值差异都敏感。总体 F 很容易捕捉本案例较大的主效应，而交互统计量只聚焦较难估计的差异之差，因此两个 P 值不必相近。模型交互 t 检验又针对模型中的平均交互系数并依赖误差模型；三者不能当成同一检验的重复验证。

10 功效与样本量：从主要 estimand 倒推设计

10.1 一元 ANOVA 近似

若仅把四个组合当作四组，可用 power.anova.test() 做初步筛查。生成机制下四个真实均值为 4.0、7.0、6.2、12.0 mmHg；下面把个体误差方差 $4.5^2$ 和诊所差异方差 $2.2^2$ 粗略相加。这个一元近似仍错误地把同诊所参与者当作独立，既不利用区组带来的精度，也不针对交互作用，因此只能用于量级检查，不能作为最终样本量方案。

true_cell_means <- c(C0_M0 = 4.0, C1_M0 = 7.0, C0_M1 = 6.2, C1_M1 = 12.0)
anova_power_at_48 <- power.anova.test(
  groups = 4, n = 48,
  between.var = var(true_cell_means), within.var = 4.5^2 + 2.2^2,
  sig.level = 0.05
)
anova_n_for_80 <- power.anova.test(
  groups = 4, power = 0.80,
  between.var = var(true_cell_means), within.var = 4.5^2 + 2.2^2,
  sig.level = 0.05
)
power_approximation <- data.frame(
  问题 = c("每组 48 人时的总体四组检验功效", "目标功效 80% 时每组所需人数"),
  估计 = c(anova_power_at_48[["power"]], ceiling(anova_n_for_80[["n"]]))
)
knitr::kable(power_approximation, digits = 3, caption = "把诊所差异粗略并入独立误差的一元 ANOVA 初步功效近似")

把诊所差异粗略并入独立误差的一元 ANOVA 初步功效近似
问题	估计
每组 48 人时的总体四组检验功效	1
目标功效 80% 时每组所需人数	10

总体四组 F 检验有功效，不代表交互作用有同样功效。样本量应围绕最重要且通常最难检测的对比计算，并预留失访、诊所退出、方差不确定性和多重主要假设的影响。

10.2 析因交互作用的模拟功效

下面在 12 家诊所、每诊所每 cell 分别 2、4、6、8 人的候选设计下，复刻诊所变异、个体误差和 RCBD 分析；每种设计模拟 400 次，以交互项双侧 $\alpha=0.05$ 为判定标准。

simulate_interaction_power <- function(per_cell, nsim = 400) {
  sim_design <- expand.grid(
    replicate = seq_len(per_cell),
    combination = seq_len(4),
    clinic_number = seq_len(n_clinics)
  )
  sim_design[["coaching"]] <- factor(c(0, 1, 0, 1)[sim_design[["combination"]]],
                                     levels = 0:1, labels = c("No", "Yes"))
  sim_design[["home_monitoring"]] <- factor(c(0, 0, 1, 1)[sim_design[["combination"]]],
                                            levels = 0:1, labels = c("No", "Yes"))
  sim_design[["clinic"]] <- factor(sim_design[["clinic_number"]])
  coaching_num <- as.integer(sim_design[["coaching"]] == "Yes")
  monitoring_num <- as.integer(sim_design[["home_monitoring"]] == "Yes")

  p_values <- numeric(nsim)
  for (s in seq_len(nsim)) {
    clinic_noise <- rnorm(n_clinics, 0, 2.2)
    y_sim <- 4 + 3 * coaching_num + 2.2 * monitoring_num +
      2.8 * coaching_num * monitoring_num +
      clinic_noise[sim_design[["clinic_number"]]] + rnorm(nrow(sim_design), 0, 4.5)
    fit_sim <- lm(y_sim ~ clinic + coaching * home_monitoring, data = sim_design)
    p_values[s] <- coef(summary(fit_sim))[
      "coachingYes:home_monitoringYes", "Pr(>|t|)"
    ]
  }
  mean(p_values < 0.05)
}

set.seed(20260916)
candidate_per_cell <- c(2, 4, 6, 8)
simulated_power <- vapply(candidate_per_cell, simulate_interaction_power, numeric(1))
power_results <- data.frame(
  每诊所每组合人数 = candidate_per_cell,
  总样本量 = n_clinics * 4 * candidate_per_cell,
  交互作用模拟功效 = simulated_power,
  Monte_Carlo标准误 = sqrt(simulated_power * (1 - simulated_power) / 400)
)
knitr::kable(power_results, digits = 3, caption = "RCBD 交互作用的模拟功效")

RCBD 交互作用的模拟功效
每诊所每组合人数	总样本量	交互作用模拟功效	Monte_Carlo标准误
2	96	0.348	0.024
4	192	0.565	0.025
6	288	0.760	0.021
8	384	0.843	0.018

set_doe_plot_font()
plot(
  power_results[["总样本量"]], power_results[["交互作用模拟功效"]],
  type = "b", pch = 19, lwd = 2, ylim = c(0, 1),
  col = palette_doe["teal"], xlab = "总样本量",
  ylab = "交互作用检验功效", main = "功效取决于最关键的对比"
)
abline(h = 0.80, lty = 2, col = palette_doe["vermillion"])
text(power_results[["总样本量"]], power_results[["交互作用模拟功效"]],
     labels = sprintf("%.2f", power_results[["交互作用模拟功效"]]), pos = 3)

折线图显示每诊所每处理组合二、四、六、八人时交互作用检验的模拟功效，并以虚线标出百分之八十目标。

不同每诊所每组合样本数下的模拟交互作用检验功效

这只是条件于效应 2.8、误差 SD 4.5、12 家诊所均留在研究中且模型正确的功效。正式设计应对更保守效应、方差估计误差、不同相关结构、缺失和诊所退出做情景分析，并报告模拟次数与 Monte Carlo 误差。

11 缺失、不依从与 ITT

11.1 分配之后发生的事情不能改写随机化

ITT 比较最初分配组，不论参与者是否完全接受干预。它保留随机化比较，通常对应“提供该策略”的政策效果。按实际接受处理做简单比较会打破随机化，因为接受程度可能由健康、动机或预后决定。

下面人为加入不依从和结局缺失，只为展示审计结构。观察概率被设置为与完整模拟结局和分组有关，因此完整案例分析可能偏离完整数据估计；真实研究中未观察结局不可用于直接建模缺失机制。

set.seed(20260917)
received_coaching_num <- rbinom(
  nrow(doe_data), 1,
  ifelse(doe_data[["coaching_num"]] == 1, 0.86, 0.08)
)
followup_probability <- plogis(
  2.2 - 0.12 * (doe_data[["bp_reduction"]] - mean(doe_data[["bp_reduction"]])) +
    0.20 * doe_data[["coaching_num"]]
)
followup_observed <- rbinom(nrow(doe_data), 1, followup_probability)
doe_data[["bp_observed"]] <- ifelse(
  followup_observed == 1, doe_data[["bp_reduction"]], NA_real_
)

adherence_table <- table(
  分配教练 = doe_data[["coaching"]],
  实际接受教练 = factor(received_coaching_num, 0:1, c("No", "Yes"))
)
knitr::kable(adherence_table, caption = "教学用分配—接受交叉表")

教学用分配—接受交叉表
	No	Yes
No	91	5
Yes	6	90

missing_table <- aggregate(
  followup_observed ~ coaching + home_monitoring, data = doe_data, FUN = mean
)
names(missing_table)[3] <- "结局观察比例"
knitr::kable(missing_table, digits = 3, caption = "按随机分组的结局观察比例")

按随机分组的结局观察比例
coaching	home_monitoring	结局观察比例
No	No	0.938
Yes	No	0.938
No	Yes	0.792
Yes	Yes	0.896

full_itt_fit <- rcbd_fit
complete_case_itt <- lm(
  bp_observed ~ clinic + coaching * home_monitoring,
  data = doe_data, na.action = na.omit
)
itt_sensitivity <- data.frame(
  分析 = c("完整模拟结局的 ITT", "仅完整案例、仍按分配组"),
  教练简单效应_无监测 = c(coef(full_itt_fit)["coachingYes"],
                            coef(complete_case_itt)["coachingYes"]),
  交互作用 = c(coef(full_itt_fit)["coachingYes:home_monitoringYes"],
               coef(complete_case_itt)["coachingYes:home_monitoringYes"]),
  分析人数 = c(nobs(full_itt_fit), nobs(complete_case_itt))
)
knitr::kable(itt_sensitivity, digits = 3, caption = "缺失后完整案例结果的教学敏感性比较")

缺失后完整案例结果的教学敏感性比较
分析	教练简单效应_无监测	交互作用	分析人数
完整模拟结局的 ITT	2.81	2.25	192
仅完整案例、仍按分配组	2.09	3.35	171

真实分析应说明缺失发生在哪个环节，按分组报告原因与时间，使用与 estimand 对齐的方法（例如基于合理 MAR 模型的多重插补/似然法），并对 MNAR 做敏感性分析。若估计依从者效果，可采用有明确定义和额外假设的 IV/CACE 等方法；不能把 per-protocol 关联直接称为随机化效果。

12 伪重复与分配单位错位

12.1 样本量不是数据行数

若教练由诊所整体实施，诊所才是随机化单位；把每名患者当成独立随机单位会造成伪重复和过窄区间。同理，一名参与者测三次血压是技术/重复测量，不是三名独立参与者。

unit_audit <- data.frame(
  场景 = c(
    "本页：个人在诊所内随机", "诊所整体切换服务", "每人测三次血压",
    "同一培养皿读取十个视野", "家庭整体接受访视"
  ),
  实验单位 = c("参与者", "诊所", "参与者", "培养皿", "家庭"),
  数据行可能对应 = c("个人", "个人", "测量时点", "视野", "家庭成员"),
  主要相关性 = c("诊所区组", "诊所内聚类", "个体内重复", "培养皿内", "家庭内")
)
knitr::kable(unit_audit, caption = "先识别分配单位，再决定有效样本量和分析层级")

先识别分配单位，再决定有效样本量和分析层级
场景	实验单位	数据行可能对应	主要相关性
本页：个人在诊所内随机	参与者	个人	诊所区组
诊所整体切换服务	诊所	个人	诊所内聚类
每人测三次血压	参与者	测量时点	个体内重复
同一培养皿读取十个视野	培养皿	视野	培养皿内
家庭整体接受访视	家庭	家庭成员	家庭内

分析模型应反映随机化单位、重复测量和处理实施层级。只加一个“聚类稳健标准误”有时能修正方差，却不能修复只有两个集群、干预与集群完全混杂或处理污染等设计缺陷。

13 何时使用更复杂的实验设计

设计	适合问题	随机化/分析关键	主要边界
Cluster randomized	干预只能按学校、社区、诊所实施	以 cluster 随机；功效含 ICC；分析保留 cluster	cluster 数太少时推断脆弱
Split-plot	一个因素只能施加到较大单位，另一个可在子单位随机	两层随机化、两种误差层级	不能用单一独立误差 ANOVA
Crossover	个体可依次接受多个可逆处理	随机序列，考虑 period、carryover、washout	不适合永久效果或不稳定疾病
Latin square	需同时控制两个正交干扰来源	每处理在每行每列一次	对缺失敏感，通常假定无高阶交互
Repeated measures	同一个体在多个时点观察响应轨迹	预定时间、协方差/混合模型、个体内相关	时点不是独立重复，缺失机制重要
Stepped wedge	资源限制下 cluster 分期切换到干预	切换顺序随机；模型控制日历时间	时间趋势与干预易混淆

设计优先于补救模型如果研究尚未开始，优先增加独立 cluster 数、改善分配隐藏、安排基线测量、减少污染和标准化结局。事后再复杂的混合模型也无法创造不存在的独立重复或恢复被预知的分配。

14 非连续结局的应用

14.1 二元结局与计数结局

同一析因设计可以有不同响应分布，但分配单位与 estimand 不变。二元结局可用 logistic/binomial 模型，计数率可用带暴露 offset 的 Poisson；若以诊所集群随机或有个体重复，还要加入与设计相符的 GEE、混合模型或随机化层级推断。

set.seed(20260918)

# 教学用二元结局：12 周是否达到预设临床反应。
response_probability <- plogis(
  -0.8 + 0.45 * doe_data[["coaching_num"]] + 0.30 * doe_data[["monitoring_num"]] +
    0.35 * doe_data[["coaching_num"]] * doe_data[["monitoring_num"]] +
    doe_data[["clinic_effect"]] / 8
)
clinical_response <- rbinom(nrow(doe_data), 1, response_probability)
binary_fit <- glm(
  clinical_response ~ clinic + coaching * home_monitoring,
  family = binomial(), data = doe_data
)

# 教学用计数结局：不同观察周数内的远程护理接触次数。
exposure_weeks <- runif(nrow(doe_data), 10, 12)
contact_mean <- exposure_weeks * exp(
  -2.0 + 0.18 * doe_data[["coaching_num"]] + 0.25 * doe_data[["monitoring_num"]] +
    0.12 * doe_data[["coaching_num"]] * doe_data[["monitoring_num"]] +
    doe_data[["clinic_effect"]] / 12
)
contact_count <- rpois(nrow(doe_data), contact_mean)
count_fit <- glm(
  contact_count ~ clinic + coaching * home_monitoring + offset(log(exposure_weeks)),
  family = poisson(), data = doe_data
)

outcome_model_summary <- data.frame(
  响应 = c("二元临床反应", "远程接触计数率"),
  交互作用链接尺度估计 = c(
    coef(binary_fit)["coachingYes:home_monitoringYes"],
    coef(count_fit)["coachingYes:home_monitoringYes"]
  ),
  指数化交互 = exp(c(
    coef(binary_fit)["coachingYes:home_monitoringYes"],
    coef(count_fit)["coachingYes:home_monitoringYes"]
  )),
  链接尺度 = c("log odds", "log rate")
)
knitr::kable(outcome_model_summary, digits = 3, caption = "相同设计下不同响应模型的语法示例")

相同设计下不同响应模型的语法示例
响应	交互作用链接尺度估计	指数化交互	链接尺度
二元临床反应	0.859	2.36	log odds
远程接触计数率	0.397	1.49	log rate

logistic 模型中的指数化交互是 odds 尺度的乘法交互，Poisson 中是 rate ratio 尺度的乘法交互；它们不等同于绝对概率差或绝对率差的交互。为政策沟通，应从模型预测四个 cell 的标准化绝对风险/率及其差异，而不是只报告链接尺度系数。

15 从方案到报告：可审核工作流

15.1 推荐分析路径

写清 PICO 与 estimand：目标总体、处理策略、结局、时间点、依从规则和效应尺度。
识别实验单位：处理在哪个最小单位独立随机和实施？观察单位是否重复？
选择设计：CRD、区组、分层、cluster 或多层随机化；说明限制条件和分配比例。
预设主要对比：cell mean、平均主效应、简单效应、交互或联合组比较；限制主要假设数。
计算功效：围绕主要 estimand，使用保守效应/方差和失访情景；集群设计纳入 ICC。
保护实施：分配隐藏、盲法、干预一致性、污染监测和结局测量标准化。
冻结分析方案：建库规则、排除、缺失、异常、协变量调整、多重性和敏感性分析。
先画绝对结果：按随机组报告样本流、cell mean/风险/率和不确定性。
拟合设计一致模型：保留区组、聚类、析因结构和随机化单位，检查诊断。
透明解释：报告估计和区间，不用“显著/不显著”代替效果大小、精度与适用边界。

15.2 报告结果的句式模板

在 12 家诊所中，共随机分配 192 名参与者；每家诊所的四个组合各 4 人。主要分析遵循 ITT，并用包含诊所固定区组、健康教练、居家监测及其交互的线性模型估计 12 周收缩压下降。对 12 家诊所等权平均，教练的平均主效应为 … mmHg（95% CI …），监测为 … mmHg（95% CI …），差异之差交互为 … mmHg（95% CI …）。四个组合的调整后均值为 …。诊所内 restricted randomization test 对 Fisher sharp global null 给出的 P 值为 …（差异之差统计量）和 …（3-df 总体 F）；它们不检验平均零交互的弱零假设。结果适用于所研究的诊所和实施策略；关于缺失、污染、依从和向其他诊所推广的限制为 …。

15.3 常见错误及修正

常见错误	为什么不对	修正
把数据行数当实验单位数	技术重复/聚类不提供独立随机化	从处理分配层级定义实验单位
先看基线 P 值再决定是否调整	随机化后基线差异检验不回答混杂	预设强预后协变量；报告描述性平衡
有交互仍把主项系数叫平均主效应	treatment coding 下主项是参考水平简单效应	用 cell mean 的明确线性组合
“一组显著、另一组不显著”即交互	两个 P 值的差不是差异的检验	直接估计差异之差及区间
忽略区组或跨区组置换	分析与真实分配机制不一致	模型控制区组；区组内重随机化
P>0.05 就宣称无效	可能是效果小，也可能估计不精确	报告估计、CI、最小重要差异和功效
按实际接受处理做主要比较	依从不是随机分配	主要 ITT；依从者效果用合适因果方法
删除异常值直到结果理想	引入结果驱动分析选择	核查记录，按预设规则做敏感性分析

16 知识检查与练习

快速自测

每家诊所整体被分配一种政策、每家测 100 人。实验单位有多少：1,200 还是 12？
在 coaching * home_monitoring treatment coding 中，coachingYes 代表什么？
为什么本页的置换不能把 192 个处理标签跨诊所任意重排？
交互作用估计为 2.5 mmHg、95% CI 为 $-0.4$ 到 5.4。能否写“没有交互作用”？
为什么总体四组 ANOVA 的 80% 功效不保证交互项也有 80% 功效？

查看答案

12 个诊所；个人是观察单位，诊所才是分配单位。有效信息还取决于 cluster 数、大小和 ICC。
在 home_monitoring=No 的参考水平下，教练 Yes 相对 No 的简单效应，不是自动对两个监测水平平均的主效应。
实际随机化限制每家诊所每组合 4 人；跨诊所置换包含不可能分配，不能正确复刻随机化参考分布。
不能。数据与一定范围的负、零和正交互均相容；应报告估计和区间，并结合临床重要界值说明精度。
总体 F 检验可由任一 cell 差异驱动；交互是特定差异之差，通常标准误更大，应单独规划功效。

练习：修改并审计设计

用 adjusted_cells 计算“有监测时教练简单效应”和“无监测时教练简单效应”，验证两者之差等于交互。
把每诊所每组合人数改为 3，重跑功效模拟并给出 Monte Carlo 标准误。结论对效应设定有多敏感？
假设教练只能由诊所整体实施，而居家监测仍在个人层面随机。画出 split-plot 随机化层级，并指出两个因素使用的误差层级。
假设 20% 参与者缺失结局。分别写出在 MAR 和 MNAR 下需要的主要分析/敏感性分析，并说明为什么单纯均值填补不合适。
将结局改为“12 周达到血压控制”的二元变量。写出四个标准化风险和风险差交互，而不仅是 logistic 系数。

练习讨论要点

权重分别为 (0,0,-1,1) 与 (-1,1,0,0)；前者减后者得到 (1,-1,-1,1)，即差异之差。
在候选向量加入 3；功效是模拟拒绝比例，Monte Carlo SE 为 $\sqrt{\hat p(1-\hat p)/B}$ 。还应改变交互效应和误差 SD 做情景分析。
教练在诊所层面随机，以诊所间变异检验；监测及交互的可估计性取决于诊所内随机化和正确的两层误差结构。需要足够多诊所。
MAR 下可用含预后/缺失预测因子的多重插补或似然模型；MNAR 需 pattern-mixture、selection 或 tipping-point 等敏感性分析。单纯均值填补压低方差且忽略不确定性。
对四个组合从 binomial 模型预测绝对风险，再计算 $[p_{11}-p_{01}]-[p_{10}-p_{00}]$ ；这与 log odds 乘法交互不是同一个 estimand。

17 速查表与提交前清单

速查表

目标	关键代码/对比	解释提醒
CRD 随机化	`sample(rep(treatment, each=n))`	保存完整生成与执行审计
RCBD 分析	`lm(y ~ block + A * B)`	每区组包含所有组合
四个 cell mean	`L %*% coef(fit)`	`L` 应与目标总体加权一致
教练平均主效应	`(-.5,.5,-.5,.5)`	平均另一个因素的两水平
交互作用	`(1,-1,-1,1)`	差异之差，尺度必须说明
全部两两比较	`TukeyHSD(aov_fit)`	用于 familywise 探索，不替代预设对比
受限随机化检验	区组内重排处理	必须复刻实际允许分配
一元功效近似	`power.anova.test()`	不能代替对主要析因对比的规划
二元结局	`glm(..., family=binomial)`	报绝对风险与链接尺度效应
计数率	`glm(...+offset(log(t)), poisson)`	检查过度离散与聚类

提交前清单

## 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       cachem_1.1.0   
##  [6] knitr_1.51      htmltools_0.5.9 rmarkdown_2.31  lifecycle_1.0.5 cli_3.6.6      
## [11] sass_0.4.10     jquerylib_0.1.4 compiler_4.6.1  tools_4.6.1     evaluate_1.0.5 
## [16] bslib_0.12.0    yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0

实验设计（Design of Experiments）与应用详解

从随机化、重复与区组到析因效应、随机化推断、功效和真实世界实施

公共卫生学院学习模块

更新于 2026年08月12日

如何使用本教程

学习目标

1 从研究问题到可估计效应

1.1 DOE 的核心术语

1.2 四个因果对比

2 六个设计支柱

2.1 随机化、重复与局部控制

2.2 完全随机设计（CRD）

2.3 为什么本案例采用 RCBD

3 2×2 析因设计：一次实验回答三个问题

3.1 为什么不做两个独立实验

3.2 数据字典与生成机制

4 探索性分析：先看 cell mean，再看模型

4.1 四组样本量、均值与标准误

4.2 交互作用的几何解释

5 RCBD 线性模型与平衡设计 ANOVA

5.1 拟合预先规定的模型

5.2 平衡设计的 ANOVA 分解

6 用模型矩阵得到四个调整后 cell mean

6.1 为什么手动线性组合值得学习

6.2 预设主效应与交互对比

6.3 Tukey 全部两两比较何时适用

7 区组如何提高精度

7.1 与忽略诊所的 CRD 分析比较

8 模型诊断：设计正确仍要检查测量与模型

8.1 残差、正态性、方差与影响

9 Restricted randomization permutation test

9.1 用真实允许的分配构造零分布

10 功效与样本量：从主要 estimand 倒推设计

10.1 一元 ANOVA 近似

10.2 析因交互作用的模拟功效

11 缺失、不依从与 ITT

11.1 分配之后发生的事情不能改写随机化

12 伪重复与分配单位错位

12.1 样本量不是数据行数

13 何时使用更复杂的实验设计

14 非连续结局的应用

14.1 二元结局与计数结局

15 从方案到报告：可审核工作流

15.1 推荐分析路径

15.2 报告结果的句式模板

15.3 常见错误及修正

16 知识检查与练习

快速自测

练习：修改并审计设计

17 速查表与提交前清单

速查表

提交前清单