R语言裂区设计方差分析（代码+数据+可视化）

医学和生信笔记

发布于 2026-04-09 20:24:27

930

文章被收录于专栏：医学和生信笔记医学和生信笔记

在实际研究中，影响结果的因素往往不止一个。比如研究某种药物的疗效时，除了药物本身，患者的性别、年龄、用药剂量等都可能同时影响结果。如果我们对每个因素单独做方差分析，不仅效率低，还会遗漏一个重要信息——因素之间的交互作用。

多因素方差分析（Multi-way ANOVA）正是用来同时分析多个因素对结果变量影响的统计方法。它可以回答以下问题：

主效应：每个因素单独对结果有没有影响？
交互作用：多个因素联合起来，是否会产生”1+1≠2”的效果？

举个例子：A药和B药单独使用时效果一般，但联合使用时镇痛时间大幅延长——这就是典型的正交互作用。反之，两药合用反而效果下降，则为负交互作用。如果不考虑交互作用，单看主效应可能会得出错误的结论。

根据实验设计的不同，多因素方差分析有多种常见形式：

析因设计：所有因素的水平两两组合，全面估计主效应和交互作用
正交设计：因素和水平较多时，用正交表选取部分组合，减少实验次数
嵌套设计：某一因素的水平嵌套在另一因素之内，两者不能交叉
裂区设计：不同因素施加在不同层级的实验单位上，兼顾精度与可行性

前面已经介绍了析因设计的方差分析、正交设计的方差分析、嵌套设计的方差分析，本篇继续介绍裂区设计的方差分析。

裂区设计常见于实验条件的改变存在难度差异的场景：某些因素（主区因素）的水平切换代价大、操作困难，只能在粗粒度的实验单位（一级单位）上随机分配；另一些因素（副区因素）则可以在更细的实验单位（二级单位）内自由随机化，从而兼顾实验的可行性与精度。本篇以家兔皮肤损伤保护实验为例（全身注射药物为主区因素，局部毒素浓度为副区因素），演示了如何用aov()中的Error(id/factorB)语法为不同层级的因素指定各自的误差项，将方差分解为"一级单位间"和"二级单位间"两部分分别检验。值得注意的是，裂区设计的数据结构与两因素重复测量设计高度相似，两者的分析思路也是相通的，区分时需结合实验的随机化方式来判断。

unsetunset嵌套设计资料的方差分析unsetunset

unsetunset裂区设计资料的方差分析unsetunset

使用孙振球《医学统计学》第4版例11-7的数据。这是一个完全随机的2*2裂区设计，家兔为一级实验单位，注射部位为二级实验单位。

试验一种全身注射抗毒素对皮肤损伤的保护作用，将10只家兔随机等分两组，一组注射抗毒素，一组注射生理盐水作对照。分组后，每只家兔取甲、乙两部位，分别随机分配注射低浓度毒素和高浓度毒素，观察指标为皮肤受损直径（mm）。试做方差分析。

data11_7 <- data.frame(
  factorA = factor(rep(c("a1","a2"),each=10)),
  factorB = factor(rep(c("b1","b2"),10)),
  id = factor(rep(c(1:10),each=2)),
  y = c(15.75,19.00,15.50,20.75,15.50,18.50,17.00,20.50,16.50,20.00,
        18.25,22.25,18.50,21.50,19.75,23.50,21.50,24.75,20.75,23.75)
  )
str(data11_7)
## 'data.frame':    20 obs. of  4 variables:
##  $ factorA: Factor w/ 2 levels "a1","a2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ factorB: Factor w/ 2 levels "b1","b2": 1 2 1 2 1 2 1 2 1 2 ...
##  $ id     : Factor w/ 10 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
##  $ y      : num  15.8 19 15.5 20.8 15.5 ...
head(data11_7)
##   factorA factorB id     y
## 1      a1      b1  1 15.75
## 2      a1      b2  1 19.00
## 3      a1      b1  2 15.50
## 4      a1      b2  2 20.75
## 5      a1      b1  3 15.50
## 6      a1      b2  3 18.50

裂区设计的A因素只作用于一级实验单位，B因素只作用于二级实验单位，所以其方差分析也是由两部分组成（P183）。如果你认真观察，你会发现这这个数据结构和两因素重复测量数据结构一致。

只看数据和代码对于了解数据结构不够直观，下面画两个图，展示数据结构：

# 创建一个直观的数据布局图
library(dplyr)
library(ggplot2)

# 准备绘图数据
plot_data <- data11_7 %>%
  mutate(
    factorA_label = factor(ifelse(factorA == "a1", "全身: 抗毒素", "全身: 生理盐水")),
    factorB_label = factor(ifelse(factorB == "b1", "局部: 低浓度", "局部: 高浓度"),
                          levels = c("局部: 低浓度", "局部: 高浓度")),
    id_label = paste("家兔", id)
  )

# 创建热图风格的数据视图
ggplot(plot_data, aes(x = factorB_label, y = reorder(id_label, as.numeric(id)))) +
  geom_tile(aes(fill = y), color = "white", size = 1) +
  geom_text(aes(label = round(y, 2)), color = "black", size = 4) +
  facet_grid(. ~ factorA_label, scales = "free", space = "free") +
  scale_fill_gradient(low = "#e3f2fd", high = "#1565c0", name = "受损直径(mm)") +
  labs(
    title = "裂区设计数据结构示意图",
    subtitle = "10只家兔 × 2个部位 = 20个观测值",
    x = "局部处理（家兔内因子）",
    y = "家兔编号"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    plot.subtitle = element_text(hjust = 0.5, color = "darkgray"),
    axis.text.x = element_text(angle = 0, hjust = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    strip.background = element_rect(fill = "#f5f5f5", color = "gray"),
    strip.text = element_text(face = "bold")
  )

# 更详细的每只家兔内部结构
library(patchwork)
## 
## Attaching package: 'patchwork'
## The following object is masked from 'package:MASS':
## 
##     area

# 为每只家兔创建一个小图
plot_list <- list()

for(i in1:10) {
  rabbit_data <- data11_7 %>% filter(id == i)

  p <- ggplot(rabbit_data, aes(x = factorB, y = y, fill = factorA)) +
    geom_bar(stat = "identity", width = 0.6) +
    geom_text(aes(label = round(y, 2)), vjust = -0.5, size = 3) +
    ylim(0, max(data11_7$y) * 1.1) +
    labs(
      title = paste("家兔", i, "-", 
                   ifelse(unique(rabbit_data$factorA) == "a1", 
                          "抗毒素组", "生理盐水组")),
      x = "局部毒素浓度",
      y = "受损直径(mm)"
    ) +
    theme_minimal() +
    theme(
      plot.title = element_text(size = 9),
      axis.text = element_text(size = 8),
      legend.position = "none"
    )

  plot_list[[i]] <- p
}

# 排列所有小图
wrap_plots(plot_list, ncol = 5) +
  plot_annotation(
    title = "每只家兔的观测结构",
    subtitle = "每只家兔接受1种全身处理 + 2种局部处理"
  )

该例题中每只家兔对应着B因素（毒素浓度）的两个水平（每只家兔会注射两种浓度的毒素），但每只家兔只对应A因素的1个水平（每只家兔只会注射一种药物，不会同时注射两种药物，要么注射抗毒素，要么注射生理盐水），所以需要为B因素指定误差项。

# factorB is nested in id，每个id对应多个factorB
# factorA和factorB有交叉，但是id只和factorB有交叉
f <- aov(y ~ factorA * factorB + Error(id/factorB), data = data11_7)
summary(f)
## 
## Error: id
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## factorA    1  63.01   63.01   28.01 0.000735 ***
## Residuals  8  18.00    2.25                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: id:factorB
##                 Df Sum Sq Mean Sq F value   Pr(>F)    
## factorB          1  63.01   63.01  252.05 2.48e-07 ***
## factorA:factorB  1   0.11    0.11    0.45    0.521    
## Residuals        8   2.00    0.25                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

结果同课本相同。第一部分是A因素主效应和误差，第二部分是：B因素主效应、A和B的交互效应、误差。

# 看下每个因素下的平均受损直径
tapply(data11_7$y, list(data11_7$factorA,data11_7$factorB),mean)
##       b1    b2
## a1 16.05 19.75
## a2 19.75 23.15

在低浓度下（b1）：
- 抗毒素组（a1）平均受损直径为：16.05
- 生理盐水组（a2）平均受损直径为：19.75
在高浓度下（b2）：
- 抗毒素组（a1）平均受损直径为：19.75
- 生理盐水组（a2）平均受损直径为：23.15

结论为：无论是低浓度毒素还是高浓度毒素所致的皮肤损伤，全身注射抗毒素的皮肤受损直径（mm）均小于对照组。全身注射抗毒素对皮肤损伤有保护作用。

裂区设计和嵌套设计R方差分析实现的参考链接：

Crossed and Nested Factors
aov() error term in R: what’s the difference bw Error(id) and Error(id/timevar) specification?
Formulae in R
R and Analysis of Variance

医学和生信笔记，专注R语言在医学中的使用！

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2026-04-06，如有侵权请联系 cloudcommunity@tencent.com 删除

数据结构