R语言ANOVA and Plot 代码

该代码使用 R markdown  语法
---
title: "ANOVA and Plot"
output:
  html_document:
    df_print: paged
---

## anova oneway单因素方差分析
```{r}
# 加载dplyr包
library(dplyr) 
# 加载cholesterol数据集
data(cholesterol, package="multcomp") 
# 计算各治疗组的样本量、均值、标准差和95%置信区间
plotdata <- cholesterol %>%
  group_by(trt) %>%
  summarize(n = n(),
            mean = mean(response),
            sd = sd(response),
            ci = qt(0.975, df = n - 1) * sd / sqrt(n)) 
# 显示计算结果
plotdata
# 进行单因素方差分析
fit <- aov(response ~ trt, data=cholesterol) 
# 查看方差分析结果
summary(fit) 
# 加载ggplot2包
library(ggplot2) 
# 绘制治疗组均值和置信区间的图形
ggplot(plotdata, aes(x = trt, y = mean, group = 1)) +
  geom_point(size = 3, color="red") +
  geom_line(linetype="dashed", color="darkgrey") +
  geom_errorbar(aes(ymin = mean - ci, ymax = mean + ci), width=.1) +
  theme_bw() +
  labs(x="Treatment", y="Response", title="Mean Plot with 95% Confidence Interval")
```
## one way and two way ANOVA
```{r}
# 加载必要的包
library(dplyr)
library(ggplot2)
library(car)
library(multcomp)

# 单因素方差分析绘图示例 - 使用cholesterol数据集
data(cholesterol, package = "multcomp")

# 1. 绘制均值图和置信区间
plotdata <- cholesterol %>%
  group_by(trt) %>%
  summarize(
    n = n(),
    mean = mean(response),
    sd = sd(response),
    se = sd / sqrt(n),
    ci = qt(0.975, df = n - 1) * se
  )

# 基础均值图
p1 <- ggplot(plotdata, aes(x = trt, y = mean, group = 1)) +
  geom_point(size = 3, color = "red") +
  geom_line(linetype = "dashed", color = "darkgrey") +
  geom_errorbar(aes(ymin = mean - ci, ymax = mean + ci), width = 0.1) +
  theme_bw() +
  labs(
    x = "Treatment",
    y = "Response",
    title = "Mean Plot with 95% Confidence Intervals"
  )

# 2. 绘制箱线图展示组间分布
p2 <- ggplot(cholesterol, aes(x = trt, y = response)) +
  geom_boxplot(fill = "lightblue", color = "black") +
  theme_bw() +
  labs(
    x = "Treatment",
    y = "Response",
    title = "Boxplot of Treatment Groups"
  )

# 3. 绘制散点图加抖动以显示原始数据点
p3 <- ggplot(cholesterol, aes(x = trt, y = response)) +
  geom_jitter(width = 0.2, alpha = 0.5, color = "blue") +
  theme_bw() +
  labs(
    x = "Treatment",
    y = "Response",
    title = "Scatter Plot with Jitter"
  )

# 4. 两因素方差分析交互效应图示例（修正版）
set.seed(123)  # 设置随机种子以确保结果可重现

# 创建因子组合（6组）
mydata <- expand.grid(
  A = factor(c("Low", "High")),
  B = factor(c("Control", "Treatment1", "Treatment2"))
)

# 为每组定义均值（设置A和B的主效应以及交互效应）
group_means <- matrix(
  c(20, 25, 30,  # A=Low时B各水平的均值
    35, 45, 50), # A=High时B各水平的均值（注意Treatment1的均值增加更多，显示交互效应）
  nrow = 2, byrow = TRUE,
  dimnames = list(c("Low", "High"), c("Control", "Treatment1", "Treatment2"))
)

# 将矩阵转换为长格式的均值向量
mydata$mean <- as.vector(group_means)

# 为每个组生成10个观测值（共60个）
mydata <- mydata[rep(1:nrow(mydata), each = 10), ]  # 复制行形成60行
mydata$y <- mydata$mean + rnorm(nrow(mydata), mean = 0, sd = 3)  # 添加随机噪声

# 拟合两因素方差分析模型
model <- aov(y ~ A * B, data = mydata)

# 使用effect包计算预测值
if (!require(effects)) install.packages("effects")
library(effects)

# 计算A和B的交互效应
effect_data <- as.data.frame(effect("A:B", model))

# 绘制交互效应图
p4 <- ggplot(effect_data, aes(x = B, y = fit, group = A, color = A)) +
  geom_line(size = 1.2) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1, size = 0.8) +
  theme_bw() +
  labs(
    x = "Factor B",
    y = "Predicted Values",
    title = "Interaction Plot: A x B",
    color = "Factor A"
  ) +
  scale_color_manual(values = c("Low" = "blue", "High" = "red"))

# 5. 残差诊断图 - 检查方差分析假设
fit <- aov(response ~ trt, data = cholesterol)
par(mfrow = c(2, 2))
plot(fit)
par(mfrow = c(1, 1))

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