R语言析因设计方差分析及主效应和单独效应计算

医学和生信笔记

发布于 2026-03-17 16:13:07

830

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

关注公众号，发送R语言或python，可获取资料

💡专注R语言在🩺生物医学中的使用

介绍如何使用R语言进行析因设计方差分析以及如何使用emmeans计算主效应。

本文目录

方差分析
主效应计算
可视化数据

unsetunset方差分析unsetunset

将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种类型的平方和也是一样的
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的结果会给出截距项，在方差分析时可忽略）。分别给出了A因素、B因素、AB交互作用、个体间的自由度、离均差平方和、均方误差、F值、P值，可以看到结果和书中（表11-5）是一致的。

说明 aov(y ~ x1 * x2, data = data11_1)等价于aov(y ~ x1 + x2 + x1:x2, data = data11_1)，表示x1的主效应、x2的主效应、x1和x2的交互效应。

car::Anova()用法：

# 使用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因素主效应达到0.01<P<0.05，拒绝H，接受H。结合样本均数的比较结果，A因素的主效应为6%，AB的交互作用为2%，均不具有统计学意义。仅B因素（缝合后时间）的主效应22%有统计学意义。

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

unsetunset主效应计算unsetunset

如何计算各因素的主效应呢？可以借助emmeans包。

library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
# x1的主效应
emmeans(f1, "x1")
## NOTE: Results may be misleading due to involvement in interactions
##      x1   emmean   SE df lower.CL upper.CL
##  束膜缝合     40 5.48 16     28.4     51.6
##  外膜缝合     34 5.48 16     22.4     45.6
## 
## Results are averaged over the levels of: x2 
## Confidence level used: 0.95

# x2的主效应
emmeans(f1, "x2")
## NOTE: Results may be misleading due to involvement in interactions
##      x2    emmean   SE df lower.CL upper.CL
##  缝合1个月     26 5.48 16     14.4     37.6
##  缝合2个月     48 5.48 16     36.4     59.6
## 
## Results are averaged over the levels of: x1 
## Confidence level used: 0.95

x1主效应为：40-34=6，解释为：在不考虑缝合时间的影响下，外膜缝合比束膜缝合的轴突通过率提高了6%；
x2主效应为：48-26=22，解释为：在不考虑缝合方法的影响下，缝合后2个月比缝合后1个月的轴突通过率提高了22%。

本例两个因素可认为不存在交互作用，若存在交互作用，主效应可能会不准（不能反映真实情况），此时应分析各因素的单独效应。

如何计算各个单元格的均值？

# 计算各个单元格的均值
emmeans(f1, c("x1","x2"))
##      x1       x2    emmean   SE df lower.CL upper.CL
##  束膜缝合 缝合1个月     28 7.75 16    11.58     44.4
##  外膜缝合 缝合1个月     24 7.75 16     7.58     40.4
##  束膜缝合 缝合2个月     52 7.75 16    35.58     68.4
##  外膜缝合 缝合2个月     44 7.75 16    27.58     60.4
## 
## Confidence level used: 0.95

缝合时间为1个月时，x1（缝合方法）的单独效应为：28-24=4；
缝合时间为2个月时，x1（缝合方法）的单独效应为：52-44=8；
束膜缝合时，x2（缝合时间）的单独效应为：52-28=24；
外膜缝合时，x2（缝合时间）的单独效应为：44-24=20。

这个结果可以直接保存为数据框，然后画图：

plot_df <- as.data.frame(emmeans(f1, c("x1","x2")))
plot_df
##      x1       x2    emmean       SE df lower.CL upper.CL
##  束膜缝合 缝合1个月     28 7.745967 16 11.57928 44.42072
##  外膜缝合 缝合1个月     24 7.745967 16  7.57928 40.42072
##  束膜缝合 缝合2个月     52 7.745967 16 35.57928 68.42072
##  外膜缝合 缝合2个月     44 7.745967 16 27.57928 60.42072
## 
## Confidence level used: 0.95

library(ggplot2)
ggplot(plot_df,aes(x1,emmean))+
  geom_line(aes(color=x2,group=x2),linewidth=2)+
  geom_point(aes(color=x2),size=6)+
  theme_bw()

最好是在一开始就把分类变量变成factor，然后规定好level，这样画图也方便。

unsetunset可视化数据unsetunset

简单介绍一下可视化两因素析因设计的其他方法：

# R自带函数，无需加载R包
interaction.plot(data11_1$x2, data11_1$x1, data11_1$y, type = "b", 
                 col = c("red","blue"), pch = c(12,15),
                 xlab = "缝合时间", ylab = "轴突通过率")

另外一种可视化方法：

library(gplots)
attach(data11_1)
plotmeans(y ~ interaction(x1,x2),
          connect = list(c(1,3), c(2,4)),
          col = c("red","darkgreen"),
          main = "两因素析因设计",
          xlab = "时间和方法的交互")

再介绍一种方法：

library(HH)

interaction2wt(y ~ x1 * x2)

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

关注我，不迷路：

三连一下，感谢支持

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

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

笔记

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

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

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

登录后参与评论

0 条评论

热度