一文搞定！R语言多因素方差分析（代码+可视化+结果解读）

医学和生信笔记

发布于 2026-04-16 15:30:53

460

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

本文目录：

2x2两因素析因设计资料的方差分析
IxJ两因素析因设计资料的方差分析
IxJxK三因素析因设计资料的方差分析
正交设计资料的方差分析
嵌套设计资料的方差分析
裂区设计资料的方差分析

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

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

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

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

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

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

之前给大家详细介绍了如何使用R语言自带的aov()函数做多因素方差分析：

今天给给大家介绍如何使用car::Anova()做多因素方差分析。关于各种研究设计类型的简介，请参考上述推文。

unsetunset2x2两因素析因设计资料的方差分析unsetunset

使用孙振球《医学统计学》第4版例11-1的数据。

将20只家兔随机等分4组，每组5只，进行神经损伤后的缝合试验。处理由两个因素组合而成，A因素为缝合方法，有两个水平，一个水平为外膜缝合，记作a1，另一个水平为束膜缝合，记作a2。B因素为缝合后的时间，有两个水平，一个水平为缝合后1个月，记作b1，另一个水平为缝合后2个月，记作b2。试验结果为每只家免神经缝合后的轴突通过率（%）。欲比较不同缝合方法及缝合后时间对轴突通过率的影响，试做析因设计的方差分析。

data11_1 <- data.frame(
  x1 = rep(c("外膜缝合","束膜缝合"), each = 10),
  x2 = rep(c("缝合1个月","缝合2个月"), each = 5),
  y = c(10,10,40,50,10,30,30,70,60,30,10,20,30,50,30,50,50,70,60,30)
)

str(data11_1)
## 'data.frame':    20 obs. of  3 variables:
##  $ x1: chr  "外膜缝合" "外膜缝合" "外膜缝合" "外膜缝合" ...
##  $ x2: chr  "缝合1个月" "缝合1个月" "缝合1个月" "缝合1个月" ...
##  $ y : num  10 10 40 50 10 30 30 70 60 30 ...
table(data11_1$x1,data11_1$x2)
##           
##            缝合1个月 缝合2个月
##   束膜缝合         5         5
##   外膜缝合         5         5

数据一共3列，第1列是缝合方法，第2列是时间，第3列是轴突通过率。

进行析因设计资料的方差分析（考虑所有因素的主效应和交互作用）：

# 使用3型平方和时一定要加上这句代码，不然可能得到错误的结果！
options(contrasts = c("contr.sum", "contr.poly"))
library(car)
## Loading required package: carData
fit <- lm(y ~ x1 * x2, data = data11_1) 
car::Anova(fit,type = 3) #计算3型平方和，变量顺序和结果无关
## Anova Table (Type III tests)
## 
## Response: y
##             Sum Sq Df F value    Pr(>F)    
## (Intercept)  27380  1 91.2667 5.167e-08 ***
## x1             180  1  0.6000   0.44987    
## x2            2420  1  8.0667   0.01182 *  
## x1:x2           20  1  0.0667   0.79955    
## Residuals     4800 16                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A因素主效应所对应的检验假设为H0：A因素主效应=0；B因素主效应所对应的检验假设为H0：B因素主效应=0；AB交互作用所对应的检验假设为H0：AB交互作用=0。

只有B因素主效应的P值在0.01到0.05之间（P<0.05），仅B因素（缝合后时间）的主效应有统计学意义。

结论为：尚不能认为两种缝合方法对神经轴突通过率有影响；可以认为缝合后2个月与1个月相比，神经轴突通过率提高了。

本例中AB交互作用的P值为0.7995，远大于0.05，因此可以认为两个因素之间不存在交互作用，缝合方法的效果不会因缝合时间不同而改变，此时直接报告和解释主效应是合适的。

💡如果交互作用显著，该怎么办？

当交互作用存在时，主效应（边际均值之差）会把不同条件下的效果”平均掉”，可能掩盖甚至歪曲真实情况，此时不建议单独解读主效应。正确的做法是转而分析简单效应（Simple Effect，也叫单独效应）——即固定某一因素在某个水平上，再看另一个因素的效果。

关于主效应、交互效应、简单效应的解释以及后续分析（简单效应分析），请参考R语言析因设计方差分析，本篇不再赘述。

unsetunsetIxJ两因素析因设计资料的方差分析unsetunset

使用孙振球《医学统计学》第4版例11-2的数据。

观察A、B两种镇痛药物联合运用在产妇分娩时的镇痛效果。A药取3个剂量：1.0mg、2.5mg、5.0mg；B药也取3个剂量：5μg、15μg、30μg。共9个处理组。将27名产妇随机等分为9组，每组3名产妇，记录每名产妇分娩时的镇痛时间。试分析A、B两药联合运用的镇痛效果。

data11_2 <- data.frame(
  druga = rep(c("1mg","2.5mg","5mg"), each = 3),
  drugb = rep(c("5微克","15微克","30微克"),each = 9),
  y = c(105,80,65,75,115,80,85,120,125,115,105,80,125,130,90,65,
        120,100,75,95,85,135,120,150,180,190,160)
)

str(data11_2)
## 'data.frame':    27 obs. of  3 variables:
##  $ druga: chr  "1mg" "1mg" "1mg" "2.5mg" ...
##  $ drugb: chr  "5微克" "5微克" "5微克" "5微克" ...
##  $ y    : num  105 80 65 75 115 80 85 120 125 115 ...
head(data11_2)
##   druga drugb   y
## 1   1mg 5微克 105
## 2   1mg 5微克  80
## 3   1mg 5微克  65
## 4 2.5mg 5微克  75
## 5 2.5mg 5微克 115
## 6 2.5mg 5微克  80
table(data11_2$druga,data11_2$drugb)
##        
##         15微克 30微克 5微克
##   1mg        3      3     3
##   2.5mg      3      3     3
##   5mg        3      3     3

数据一共3列，第1列是a药物的剂量（3种剂量，代表3个水平），第2列是b药物的剂量（3种剂量），第3列是镇痛时间。

进行两因素三水平的析因设计资料方差分析（考虑所有因素的主效应和交互作用）：

# 使用3型平方和时一定要加上这句代码，不然可能得到错误的结果！
options(contrasts = c("contr.sum", "contr.poly"))
library(car)
fit <- lm(y ~ druga * drugb, data = data11_2) 
car::Anova(fit,type = 3)#计算3型平方和，变量顺序和结果无关
## Anova Table (Type III tests)
## 
## Response: y
##             Sum Sq Df  F value    Pr(>F)    
## (Intercept) 326700  1 842.0907 < 2.2e-16 ***
## druga         6572  2   8.4702  0.002556 ** 
## drugb         7022  2   9.0501  0.001905 ** 
## druga:drugb   7872  4   5.0728  0.006467 ** 
## Residuals     6983 18                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

结果和表11-9也是一模一样的！

结论：①A药不同剂量的镇痛效果不同；②B药不同剂量的镇痛效果不同；③A、B两药有交互作用，A药5.0mg和B药30μg时，镇痛时间持续最长。

unsetunsetIxJxK三因素析因设计资料的方差分析unsetunset

使用孙振球《医学统计学》第4版例11-3的数据。

用5×2×2析因设计研究5种类型的军装在两种环境、两种活动状态下的散热效果，将100名受试者随机等分20组，观察指标是受试者的主观热感觉（从“冷”到“热”按等级评分）。试进行方差分析。

data11_3 <- foreign::read.spss("../datasets/例11-03-5种军装热感觉5-2-2.sav", 
                             to.data.frame = T)

data11_3$a <- factor(data11_3$a)

str(data11_3)
## 'data.frame':    100 obs. of  4 variables:
##  $ b: Factor w/ 2 levels "干燥","潮湿": 1 1 1 1 1 1 1 1 1 1 ...
##  $ c: Factor w/ 2 levels "静坐","活动": 1 1 1 1 1 1 1 1 1 1 ...
##  $ a: Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 2 2 2 2 2 ...
##  $ x: num  0.25 -0.25 1.25 -0.75 0.4 ...
##  - attr(*, "variable.labels")= Named chr [1:4] "活动环境" "活动状态" "军装类型" "主观热感觉"
##   ..- attr(*, "names")= chr [1:4] "b" "c" "a" "x"
##  - attr(*, "codepage")= int 65001
head(data11_3)
##      b    c a     x
## 1 干燥 静坐 1  0.25
## 2 干燥 静坐 1 -0.25
## 3 干燥 静坐 1  1.25
## 4 干燥 静坐 1 -0.75
## 5 干燥 静坐 1  0.40
## 6 干燥 静坐 2  0.30

a：军装类型，5个类型
b：活动环境：干燥和潮湿
c：活动状态：静坐和活动

进行3因素析因设计资料的方差分析（考虑所有的主效应和交互作用）：

# 使用3型平方和时一定要加上这句代码，不然可能得到错误的结果！
options(contrasts = c("contr.sum", "contr.poly"))
library(car)
fit <- lm(x ~ a * b * c, data = data11_3) 
car::Anova(fit,type = 3)#计算3型平方和，变量顺序和结果无关
## Anova Table (Type III tests)
## 
## Response: x
##             Sum Sq Df   F value    Pr(>F)    
## (Intercept) 660.85  1 1538.1013 < 2.2e-16 ***
## a             5.20  4    3.0245   0.02238 *  
## b             9.94  1   23.1382 6.977e-06 ***
## c           283.35  1  659.4854 < 2.2e-16 ***
## a:b           1.94  4    1.1284   0.34909    
## a:c           1.48  4    0.8620   0.49053    
## b:c          12.68  1   29.5139 5.824e-07 ***
## a:b:c         1.61  4    0.9366   0.44718    
## Residuals    34.37 80                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

结果也是和表11-12一模一样。

结论：不同军装(a)、不同环境(b)和不同活动状态(c)的主观热感觉的主效应都有差别，但尚不能认为军装类型的主观热感觉与其他两个试验因素（环境、活动状态）存在交互作用(a:b,a:c,a：b:c)。

进一步计算不同军装类型(a)的主观热感觉小计（加和）可得：第4种类型的军装具有散热效果，第5种类型的军装具有保温效果。

# 不同军装类型的主观热感觉小计
tapply(data11_3$x,data11_3$a,sum)
##     1     2     3     4     5 
## 51.80 52.22 51.10 43.81 58.14

unsetunset正交设计资料的方差分析unsetunset

正交设计是析因设计的一种缩减方案——在因素和水平较多时，全面析因设计的实验次数会急剧增加，正交设计通过选取一部分有代表性的组合（正交表），在保证主效应可估的前提下大幅减少实验次数。正交设计的核心特点是用部分实验组合（正交表）来估计主效应，同时牺牲高阶交互作用的信息。

使用孙振球《医学统计学》第4版例11-4的数据。

研究雌螺产卵的最优条件，在20cm²的泥盒里饲养同龄雌螺10只，试验条件有4个因素，每个因素2个水平。试在考虑温度与含氧量对雌螺产卵有交互作用的情况下安排正交试验。

data11_4 <- data.frame(
  a = rep(c("5度","25度"),each = 4),
  b = rep(c(0.5, 5.0), each = 2),
  c = c(10, 30),
  d = c(6.0, 8.0,8.0,6.0,8.0,6.0,6.0,8.0),
  x = c(86,95,91,94,91,96,83,88)
)

data11_4$a <- factor(data11_4$a)
data11_4$b <- factor(data11_4$b)
data11_4$c <- factor(data11_4$c)
data11_4$d <- factor(data11_4$d)

str(data11_4)
## 'data.frame':    8 obs. of  5 variables:
##  $ a: Factor w/ 2 levels "25度","5度": 2 2 2 2 1 1 1 1
##  $ b: Factor w/ 2 levels "0.5","5": 1 1 2 2 1 1 2 2
##  $ c: Factor w/ 2 levels "10","30": 1 2 1 2 1 2 1 2
##  $ d: Factor w/ 2 levels "6","8": 1 2 2 1 2 1 1 2
##  $ x: num  86 95 91 94 91 96 83 88
head(data11_4)
##      a   b  c d  x
## 1  5度 0.5 10 6 86
## 2  5度 0.5 30 8 95
## 3  5度   5 10 8 91
## 4  5度   5 30 6 94
## 5 25度 0.5 10 8 91
## 6 25度 0.5 30 6 96

进行正交设计资料的方差分析，只考虑4个因素的主效应以及a和b的一阶交互作用：

# 使用3型平方和时一定要加上这句代码，不然可能得到错误的结果！
options(contrasts = c("contr.sum", "contr.poly"))
library(car)
fit <- lm(x ~ a + b + c + d + a:b, data = data11_4) 
car::Anova(fit,type = 3)#计算3型平方和，变量顺序和结果无关
## Anova Table (Type III tests)
## 
## Response: x
##             Sum Sq Df F value    Pr(>F)    
## (Intercept)  65522  1 26208.8 3.815e-05 ***
## a                8  1     3.2   0.21554    
## b               18  1     7.2   0.11535    
## c               60  1    24.2   0.03893 *  
## d                4  1     1.8   0.31175    
## a:b             50  1    20.0   0.04654 *  
## Residuals        5  2                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

结果和表11-23一模一样，结论：雌螺产卵条件主要与泥土含水量(c)、温度与含氧量的交互作用(a:b)有关。

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

嵌套设计（也称系统分组设计）用于描述因素之间存在层级包含关系的情形——二级因素的水平并不是在所有一级因素水平下都相同，而是每个一级水平下各自拥有一套独立的二级水平。换言之，两个因素之间只有”包含”关系，没有”交叉”关系，因此不能估计交互作用，只能分别估计各层级因素的主效应。

嵌套设计的方差分析也是要注意指定主效应和交互效应。

使用孙振球《医学统计学》第4版例11-6的数据。

试验甲、乙、丙3种催化剂在不同温度下对某化合物的转化作用。由于各催化剂所要求的温度范围不同，将催化剂作为一级实验因素（I=3），温度作为二级实验因素（J=3），采用嵌套设计，每个处理重复2次（n=2）。试做方差分析。

data11_6 <- data.frame(
  factor1 = factor(rep(c("A","B","C"),each=6)),
  factor2 = factor(rep(c(70,80,90,55,65,75,90,95,100),each=2)),
  y = c(82,84,91,88,85,83,65,61,62,59,56,60,71,67,75,78,85,89)
  )
str(data11_6)
## 'data.frame':    18 obs. of  3 variables:
##  $ factor1: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 2 2 2 2 ...
##  $ factor2: Factor w/ 8 levels "55","65","70",..: 3 3 5 5 6 6 1 1 2 2 ...
##  $ y      : num  82 84 91 88 85 83 65 61 62 59 ...
head(data11_6)
##   factor1 factor2  y
## 1       A      70 82
## 2       A      70 84
## 3       A      80 91
## 4       A      80 88
## 5       A      90 85
## 6       A      90 83

factor1是一级实验因素（不同的催化剂），factor2是二级实验因素（不同的温度），y是因变量。

对这个数据做个简单的可视化，方便查看研究设计结构：

library(ggplot2)

ggplot(data11_6, aes(x = factor2, y = y)) +
  # 绘制原始数据点（轻微抖动避免重叠），按催化剂分色
  geom_jitter(aes(color = factor1), width = 0.1, size = 3) +
  facet_wrap(~factor1, labeller = label_both) +
  labs(title = "嵌套设计数据可视化 (催化剂与温度)",
       x = "温度",
       y = "转化率 (%)",
       color = "催化剂") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5),
        strip.background = element_rect(fill = "lightblue"))

对于该示例这种嵌套设计，也能用car::Anova做方差分析，但是并不建议，因为无法计算3型平方和，且2型平方和的计算也存在缺陷。

# 强行使用也能得出一样的结果，但是会给出警告
# 如果选择3型平方和，则会直接报错
# 使用3型平方和时一定要加上这句代码，不然可能得到错误的结果！
#options(contrasts = c("contr.sum", "contr.poly"))
library(car)
fit <- lm(y ~ factor1 + factor1:factor2, data = data11_6) 
car::Anova(fit,type = 2)
## Note: model has aliased coefficients
##       sums of squares computed by model comparison
## Anova Table (Type II tests)
## 
## Response: y
##                 Sum Sq Df F value    Pr(>F)    
## factor1         1956.0  2 177.818 5.831e-08 ***
## factor1:factor2  401.0  6  12.152 0.0007156 ***
## Residuals         49.5  9                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

结果和表11-26相同。

结论：催化剂(factor1)影响该化合物的转化率。对于同一催化剂，不同温度(factor2)对转化率亦有影响。

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

裂区设计常见于实验条件的改变存在难度差异的场景：某些因素（主区因素）的水平切换代价大、操作困难，只能在粗粒度的实验单位（一级单位）上随机分配；另一些因素（副区因素）则可以在更细的实验单位（二级单位）内自由随机化，从而兼顾实验的可行性与精度。

本篇以家兔皮肤损伤保护实验为例（全身注射药物为主区因素，局部毒素浓度为副区因素），演示了如何用aov()中的Error(id/factorB)语法为不同层级的因素指定各自的误差项，将方差分解为”一级单位间”和”二级单位间”两部分分别检验。

值得注意的是，裂区设计的数据结构与两因素重复测量设计高度相似，两者的分析思路也是相通的，区分时需结合实验的随机化方式来判断。

使用孙振球《医学统计学》第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)

# 为每只家兔创建一个小图
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种局部处理"
  )

由于该数据和两因素两水平重复测量设计的结构完全一致：

A因素可以看作是组别因素
B因素可以看作是时间因素

裂区设计中B因素的不同水平施加在同一实验单位的不同部位（空间重复），而重复测量是同一实验单位在不同时间点测量，两者误差结构相同，所以可以用相同的分析方法，但实验背景不同。

因此对于这个数据，可以使用car::Anova做方差分析，我们先把这个数据变成宽数据：

library(tidyr)
# 先变成宽数据
data_wide <- data11_7 %>%
  select(id, factorA, factorB, y) %>%
  pivot_wider(names_from = factorB, values_from = y) %>%
  rename(b1 = `b1`, b2 = `b2`)
head(data_wide)
## # A tibble: 6 × 4
##   id    factorA    b1    b2
##   <fct> <fct>   <dbl> <dbl>
## 1 1     a1       15.8  19  
## 2 2     a1       15.5  20.8
## 3 3     a1       15.5  18.5
## 4 4     a1       17    20.5
## 5 5     a1       16.5  20  
## 6 6     a2       18.2  22.2

然后使用car::Anova进行方差分析：

# 建立多元线性模型
mlm_model <- lm(cbind(b1,b2) ~ factorA, data = data_wide)
# 定义组内因子
time_info <- data.frame(factorB=factor(c("b1","b2"),
                                       levels = c("b1","b2")))
# 使用3型平方和时一定要加上这句代码，不然可能得到错误的结果！
options(contrasts = c("contr.sum", "contr.poly"))
suppressMessages(library(car))
# 使用Anova()进行重复测量方差分析，并启用球形检验
res <- Anova(mlm_model, idata = time_info, 
             idesign = ~factorB, type = "III")
# 查看方差分析结果
summary(res, multivariate = FALSE)
## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##                 Sum Sq num Df Error SS den Df  F value    Pr(>F)    
## (Intercept)     7742.1      1       18      8 3440.939 7.923e-12 ***
## factorA           63.0      1       18      8   28.006 0.0007354 ***
## factorB           63.0      1        2      8  252.050 2.480e-07 ***
## factorA:factorB    0.1      1        2      8    0.450 0.5212275    
## ---
## 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）均小于对照组。全身注射抗毒素对皮肤损伤有保护作用。

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

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

数据结构

本文分享自医学和生信笔记微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

登录后参与评论

0 条评论

热度