R语言方差分析作图

# 加载必要的包
library(dplyr)
library(ggplot2)
library(car)
library(effects)
library(multcomp)
library(gridExtra)

# 设置随机种子以确保结果可重现
set.seed(123)

# ----------------------------
# 单因素方差分析绘图示例
# ----------------------------

# 使用内置的PlantGrowth数据集
data(PlantGrowth)

# 1. 箱线图：展示各组分布
p1 <- ggplot(PlantGrowth, aes(x = group, y = weight)) +
  geom_boxplot(fill = "lightblue", color = "black") +
  theme_minimal() +
  labs(
    title = "箱线图：不同处理组的植物重量",
    x = "处理组",
    y = "重量"
  )

# 2. 均值图 + 置信区间
summary_data <- PlantGrowth %>%
  group_by(group) %>%
  summarize(
    mean = mean(weight),
    sd = sd(weight),
    n = n(),
    se = sd / sqrt(n),
    lower = mean - qt(0.975, df = n - 1) * se,
    upper = mean + qt(0.975, df = n - 1) * se
  )

p2 <- ggplot(summary_data, aes(x = group, y = mean, group = 1)) +
  geom_point(size = 3, color = "red") +
  geom_line(linetype = "dashed", color = "darkgrey", linewidth = 0.7) +  # 将size改为linewidth
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1, linewidth = 0.8) +  # 将size改为linewidth
  theme_minimal() +
  labs(
    title = "均值图 + 95%置信区间",
    x = "处理组",
    y = "平均重量"
  )

# 3. 散点图 + 抖动
p3 <- ggplot(PlantGrowth, aes(x = group, y = weight)) +
  geom_jitter(width = 0.2, alpha = 0.6, color = "blue") +
  stat_summary(fun = "mean", geom = "point", color = "red", size = 4) +
  theme_minimal() +
  labs(
    title = "散点图 + 均值点",
    x = "处理组",
    y = "重量"
  )

# ----------------------------
# 双因素方差分析绘图示例
# ----------------------------

# 创建示例数据集
n_per_group <- 10  # 每组样本量

# 创建因子组合
mydata <- expand.grid(
  gender = factor(c("Male", "Female")),
  treatment = factor(c("Control", "Low", "High"))
)

# 为每组定义均值（包含主效应和交互效应）
means_matrix <- matrix(
  c(50, 60, 75,  # 男性在各处理水平下的均值
    55, 75, 70), # 女性在各处理水平下的均值
  nrow = 2, byrow = TRUE,
  dimnames = list(c("Male", "Female"), c("Control", "Low", "High"))
)

# 生成数据
mydata$mean <- as.vector(means_matrix)
mydata <- mydata[rep(1:nrow(mydata), each = n_per_group), ]
mydata$score <- mydata$mean + rnorm(nrow(mydata), 0, 8)

# 4. 交互效应图
model <- aov(score ~ gender * treatment, data = mydata)
effect_data <- as.data.frame(effect("gender:treatment", model))

p4 <- ggplot(effect_data, aes(x = treatment, y = fit, group = gender, color = gender)) +
  geom_line(linewidth = 1.2) +  # 将size改为linewidth
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1, linewidth = 0.8) +  # 将size改为linewidth
  theme_minimal() +
  labs(
    title = "性别 × 处理 交互效应图",
    x = "处理水平",
    y = "预测分数",
    color = "性别"
  ) +
  scale_color_manual(values = c("Male" = "blue", "Female" = "red"))

# 5. 箱线图：双因素设计
p5 <- ggplot(mydata, aes(x = treatment, y = score, fill = gender)) +
  geom_boxplot(position = position_dodge(0.8)) +
  theme_minimal() +
  labs(
    title = "双因素箱线图",
    x = "处理水平",
    y = "分数",
    fill = "性别"
  )

# 6. 散点图 + 分组
p6 <- ggplot(mydata, aes(x = treatment, y = score, color = gender)) +
  geom_jitter(alpha = 0.6, position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8)) +
  stat_summary(fun = "mean", geom = "point", size = 4, position = position_dodge(0.8)) +
  theme_minimal() +
  labs(
    title = "双因素散点图",
    x = "处理水平",
    y = "分数",
    color = "性别"
  )

# 7. 残差诊断图
fit <- aov(weight ~ group, data = PlantGrowth)
par(mfrow = c(2, 2))
plot(fit)
par(mfrow = c(1, 1))

# 显示所有图形
grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 2)