R语言方差分析“避坑指南”

医学和生信笔记

发布于 2026-06-05 21:02:56

100

方差分析那些坑

R语言做方差分析、协方差分析有非常多的注意事项，如果你不是“R语言高手+统计学高手”，那你很容易掉到这些“坑”里。今天给大家盘一盘我遇到的那些“坑”。

先说说R自带的aov的“坑”，这个主要和方差分析3种类型平方和的计算方式有关，各个方法的特点如下：

ANOVA 类型	特点
Type I（顺序平方和）（aov() 默认）	平方和按模型中变量输入顺序依次分配；先放入的变量先“拿走”能解释的变异。
Type II（部分平方和）	每个主效应在调整了其他所有主效应后计算（但不调整交互项）。适用于无交互的加性模型。
Type III（边际平方和）	每个效应在调整了模型中所有其他效应（包括交互）后计算。常用于非均衡设计或含交互模型。

aov默认使用1型平方和，且无法更改，且帮助文档中也说aov是为均衡设计准备的。

这就导致如果是非均衡数据或者有协变量等情况下，你在公式中先写哪个变量后写哪个变量对结果影响很大，你需要非常了解你的研究设计类型和研究目的以及R语言中aov函数的注意事项，才能得到正确的结果。示例请见：R语言方差分析注意事项

关于aov，我给大家总结了以下4条注意事项（默认都没有缺失值的情况下）：

1. 均衡设计且无协变量：三种平方和结果相同，aov()可放心使用。
2. 非均衡设计但单因素（只有一个处理因素）：平方和类型无影响，aov()结果可靠。
3. 非均衡设计且多因素（有多个处理因素）：aov()默认使用I型平方和，结果受变量顺序影响；如需II型或III型，应使用car::Anova()。
4. 有协变量的设计：aov()的I型平方和受变量顺序显著影响；推荐使用car::Anova()的II型或III型平方和。

但是你以为car::Anova能完美取代aov吗？并不能，只要你用的够多，就会源源不断的碰到问题……

下面是我在一个示例中发现的新“坑”。示例如下（孙振球《医学统计学》第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列是轴突通过率。

使用aov进行析因设计资料的方差分析（考虑所有因素的主效应和交互作用）不会有任何问题：

# 完全均衡的设计，调换变量顺序无影响，3种类型的平方和也是一样的
f1 <- aov(y ~ x1 * x2, data = data11_1)
summary(f1)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## x1           1    180     180   0.600 0.4499  
## x2           1   2420    2420   8.067 0.0118 *
## x1:x2        1     20      20   0.067 0.7995  
## Residuals   16   4800     300                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

但是，如果你直接用car::Anova()会得到以下结果：

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)   3920  1 13.0667 0.002325 **
## x1              40  1  0.1333 0.719783   
## x2            1440  1  4.8000 0.043610 * 
## x1:x2           20  1  0.0667 0.799545   
## Residuals     4800 16                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

你会发现两种方法得到的结果完全不一样（平方和差很多，导致F值和P值都不一样）！aov的结果是对的，car::Anova是错的！

这是因为car::Anova的默认内部对比方式不同，3型平方和对编码方式非常敏感（关于各种编码方式的区别，请看：分类变量进行回归分析时编码方式）。

如果你要使用car::Anova计算3型平方和，那么需要加上一句options(contrasts=c("contr.sum","contr.poly"))才能得到正确的结果：

# 使用3型平方和时一定要加上这句代码，不然可能得到错误的结果！
options(contrasts = c("contr.sum", "contr.poly"))
library(car)
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

默认的编码方式是contr.treatment。

# treatment对比：基线比较
# x1束膜缝合 = 相对于外膜缝合的差异
# 当有x1*x2交互时，这个"相对于基线"的定义变得模糊

# sum对比：与总体均值比较  
# x1束膜缝合 = 与总体均值的差异
# x1外膜缝合 = 与总体均值的差异（但符号相反）
# 更符合"调整所有其他效应"的理念

所以，在使用car::Anova计算3型平方和时，需要加上options(contrasts=c("contr.sum", "contr.poly"))！

除此之外，car::Anova也不适合用在嵌套设计的方差分析，示例如下。

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

# 进行方差分析
library(car)
fit <- lm(y ~ factor1 + factor1:factor2, data = data11_6) 
car::Anova(fit,type = 2) # 2型平方和会给出note提示
## 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

结果中给出了提示：“model has aliased coefficients”（模型存在别名系数）。虽然lm()函数为了能跑通模型，自动剔除（或设为NA）了部分无法估计的参数（别名系数），使得模型勉强能拟合，但这意味着模型参数并不是全都可以独立估计的。

“sums of squares computed by model comparison”（平方和通过模型比较计算）。Type II SS通常通过对比系数矩阵直接计算，但是因为参数有缺失，软件无法直接通过公式计算，而是采用了模型比较法（Model Comparison）。

计算3型平方和则会直接报错：

Type III平方和（3型平方和）的计算逻辑是“在调整了所有其他效应后，该效应的贡献”。它要求模型必须能估计出所有可能的交互作用，并且需要设置正交对比。但是嵌套数据设计矩阵中存在大量“空单元格”，这导致模型参数无法唯一估计（出现Aliased Coefficients），R会自动将部分系数标记为NA。所以无法计算。

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

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

模型