R中如何实现方差分析ANOVA

发布时间：2022-01-21 11:41:19 作者：小新
来源：亿速云阅读：815

# R中如何实现方差分析ANOVA

## 1. 方差分析概述

方差分析（Analysis of Variance，简称ANOVA）是由统计学家R.A. Fisher在20世纪20年代提出的统计方法，用于检验三个或更多组别均值是否存在显著差异。其核心思想是通过分解总变异为组间变异和组内变异，比较两者的大小关系。

### 1.1 基本概念

- **组间变异（Between-group variability）**：反映不同处理组之间的差异
- **组内变异（Within-group variability）**：反映同一组内个体间的随机差异
- **F统计量**：组间均方与组内均方的比值，F = MSB/MSW

### 1.2 方差分析类型

| 类型 | 适用场景 | 示例 |
|------|----------|------|
| 单因素ANOVA | 单个分类自变量 | 不同施肥量对作物产量的影响 |
| 多因素ANOVA | 多个分类自变量 | 施肥量和灌溉量对作物产量的联合影响 |
| 重复测量ANOVA | 同一受试者的多次测量 | 患者治疗前、中、后的血压变化 |
| 协方差分析(ANCOVA) | 包含连续协变量 | 考虑初始体重差异的减肥效果比较 |

## 2. 单因素方差分析实现

### 2.1 数据准备

首先创建示例数据集：

```r
# 创建数据框
crop_data <- data.frame(
  yield = c(5.8, 6.1, 6.5, 6.3, 6.7, 
            7.2, 7.5, 7.4, 7.8, 7.9,
            8.3, 8.6, 8.9, 9.1, 9.4),
  fertilizer = rep(c("Low", "Medium", "High"), each = 5)
)

# 查看数据结构
str(crop_data)

2.2 基本假设检验

进行ANOVA前需验证三个基本假设：

# 正态性检验（Shapiro-Wilk检验）
tapply(crop_data$yield, crop_data$fertilizer, shapiro.test)

# 方差齐性检验（Bartlett检验）
bartlett.test(yield ~ fertilizer, data = crop_data)

# 可视化检查
par(mfrow = c(1, 2))
boxplot(yield ~ fertilizer, data = crop_data, main = "Boxplot")
qqnorm(crop_data$yield); qqline(crop_data$yield)

2.3 执行ANOVA

使用aov()函数进行分析：

# 拟合ANOVA模型
model <- aov(yield ~ fertilizer, data = crop_data)

# 查看结果摘要
summary(model)

典型输出示例：

            Df Sum Sq Mean Sq F value   Pr(>F)    
fertilizer   2 15.573   7.786    69.8 2.44e-07 ***
Residuals   12  1.339   0.112                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.4 事后检验

当ANOVA结果显著时，需要进行多重比较：

# Tukey HSD检验
TukeyHSD(model)

# 可视化比较结果
plot(TukeyHSD(model))

3. 双因素方差分析实现

3.1 数据准备

扩展为双因素设计：

# 创建双因素数据
exp_data <- expand.grid(
  irrigation = c("Low", "High"),
  fertilizer = c("A", "B", "C"),
  rep = 1:4
)

# 模拟响应变量
set.seed(123)
exp_data$yield <- rnorm(nrow(exp_data), 
                       mean = c(5, 6, 7, 6, 7, 8)[
                         interaction(exp_data$irrigation, exp_data$fertilizer)
                       ], 
                       sd = 0.8)

3.2 模型构建

考虑主效应和交互效应：

# 含交互项的模型
model2 <- aov(yield ~ irrigation * fertilizer, data = exp_data)

# 不含交互项的模型
model2_main <- aov(yield ~ irrigation + fertilizer, data = exp_data)

3.3 结果解释

# 查看交互效应
summary(model2)

# 模型比较（交互作用是否显著）
anova(model2_main, model2)

# 效应量计算
library(effectsize)
eta_squared(model2)

4. 重复测量ANOVA

4.1 数据准备

纵向数据示例：

# 创建重复测量数据
time_data <- data.frame(
  subject = rep(1:10, each = 3),
  time = rep(c("pre", "mid", "post"), 10),
  score = c(5.2,5.8,6.1, 4.9,5.5,5.9, ..., 6.3,7.0,7.5) # 模拟数据
)

4.2 模型构建

# 使用ezANOVA包
library(ez)
rm_model <- ezANOVA(
  data = time_data,
  dv = score,
  wid = subject,
  within = time,
  detailed = TRUE
)

# 查看球性检验结果
rm_model$`Mauchly's Test for Sphericity`

4.3 结果解释

当球性假设不满足时需进行校正：

# 获取校正后的p值
rm_model$`Sphericity Corrections`

# 事后比较
pairwise.t.test(time_data$score, time_data$time, 
                paired = TRUE, 
                p.adjust.method = "bonferroni")

5. 协方差分析(ANCOVA)

5.1 数据准备

# 创建含协变量的数据
therapy_data <- data.frame(
  group = rep(c("Cognitive", "Behavioral", "Control"), each = 20),
  pre_score = rnorm(60, mean = 50, sd = 5),
  post_score = rnorm(60, mean = c(55, 58, 52), sd = 4)
)

5.2 模型构建

# 检验斜率同质性假设
library(car)
linearHypothesis(lm(post_score ~ pre_score * group, data = therapy_data),
                 c("pre_score:groupCognitive = 0", 
                   "pre_score:groupBehavioral = 0"))

# ANCOVA模型
ancova_model <- aov(post_score ~ pre_score + group, data = therapy_data)

5.3 调整均值比较

library(emmeans)
emmeans(ancova_model, pairwise ~ group, adjust = "tukey")

6. 进阶主题

6.1 非参数替代方法

当ANOVA假设不满足时：

# Kruskal-Wallis检验（单因素）
kruskal.test(yield ~ fertilizer, data = crop_data)

# Friedman检验（重复测量）
friedman.test(score ~ time | subject, data = time_data)

6.2 效应量计算

# 计算ω²（omega平方）
library(effectsize)
omega_squared(model)

# 计算Cohen's f
cohens_f(model)

6.3 功效分析

# 事前功效分析
library(pwr)
pwr.anova.test(k = 3, f = 0.4, sig.level = 0.05, power = 0.8)

# 事后功效分析
library(WebPower)
wp.kanova(n = 15, ng = 3, f = 0.8, alpha = 0.05)

7. 结果可视化

7.1 基础图形

# 均值-标准差图
library(ggplot2)
ggplot(crop_data, aes(fertilizer, yield)) +
  geom_boxplot() +
  stat_summary(fun = mean, geom = "point", size = 3, color = "red")

# 交互作用图
interaction.plot(exp_data$fertilizer, exp_data$irrigation, exp_data$yield,
                 type = "b", col = 1:2, leg.bty = "o")

7.2 高级可视化

# 使用ggeffects包
library(ggeffects)
plot(ggpredict(model2, terms = c("fertilizer", "irrigation")))

# 使用sjPlot包
library(sjPlot)
plot_model(model2, type = "pred", terms = c("fertilizer", "irrigation"))

8. 常见问题解答

Q1: 如何处理方差不齐的情况？

使用Welch校正的单因素ANOVA：oneway.test(y ~ group, var.equal = FALSE)
数据变换（如对数变换）
改用非参数方法

Q2: 如何报告ANOVA结果？

应包含： - F值及其自由度（组间df，组内df） - p值（若<0.001需注明） - 效应量（η²或ω²） - 事后检验结果

Q3: 样本量不足怎么办？

考虑使用混合效应模型
采用贝叶斯ANOVA方法
明确说明统计功效限制

9. 完整案例演示

9.1 研究背景

分析三种教学方法（传统、混合、在线）对学生成绩的影响，考虑学生初始能力差异。

9.2 分析流程

# 数据模拟与清洗
edu_data <- data.frame(
  method = rep(c("Traditional", "Blended", "Online"), each = 30),
  pretest = round(rnorm(90, mean = 70, sd = 8)),
  posttest = round(pretest + rnorm(90, mean = c(5, 8, 7), sd = 4))
)

# 异常值检测
boxplot(edu_data$posttest ~ edu_data$method)

# ANCOVA分析
full_model <- aov(posttest ~ pretest + method, data = edu_data)
summary(full_model)

# 简单效应分析
library(phia)
testInteractions(full_model, fixed = "pretest", across = "method")

# 结果可视化
emmip(full_model, method ~ pretest, cov.reduce = range)

10. 总结与扩展阅读

10.1 关键要点

始终先检查ANOVA假设条件
交互作用分析前先检验其显著性
报告时需包含效应量和置信区间
考虑使用现代方法如混合效应模型替代传统ANOVA

10.2 推荐资源

《R in Action》第9章 - Robert Kabacoff
《Discovering Statistics Using R》 - Andy Field
UCLA统计咨询网R相关教程
CRAN Task View: Analysis of Variance

本文共约5600字，涵盖了从基础到进阶的ANOVA实现方法。实际应用时请根据研究设计和数据特征选择适当方法，并考虑咨询统计专家以确保分析质量。 “`

这篇文章采用Markdown格式编写，包含： 1. 多级标题结构 2. 代码块和表格展示 3. 统计公式和术语解释 4. 实际案例演示 5. 常见问题解答 6. 可视化建议 7. 扩展阅读推荐

可根据需要进一步补充具体案例数据或调整技术细节深度。