lvzhaopeng的数字花园

R语言那些事 统计系列（三） rstatix and ggpubr 可视化之分面图

7分钟阅读

缘由

统计是生信的基石之一

碎碎念

要爱一个具体的人，不要爱一个抽象的人； 要爱生活，不要爱生活的意义。 ——《卡拉马佐夫兄弟》

Ref

How to Add P-values to GGPLOT Facets: Step by Step Guide - Datanovia

多分组统计分析与可视化完整指南

1. 数据加载与预处理

加载必要的R包

library(ggpubr)    # 数据可视化
library(rstatix)   # 统计检验
library(tidyverse) # 数据处理

数据准备与探索

df <- ToothGrowth
df$dose <- as.factor(df$dose)
# Add a random grouping variable
df$group <- factor(rep(c("grp1", "grp2"), 30))
head(df, 3)

数据说明：

  • ToothGrowth: R内置数据集，研究维生素C对牙齿生长的影响
  • 变量转换: 将dose从数值型转为因子型，便于分组分析
  • 新增分组: 添加随机分组变量，模拟更复杂的实验设计
  • 数据结构: 包含len(牙齿长度)、supp(补充类型)、dose(剂量)、group(分组)

2. 单分组变量的统计分析

基础统计检验

stat.test <- df %>%
	group_by(dose) %>%           # 按剂量分组
	t_test(len ~ supp) %>%       # 比较不同补充类型的牙齿长度
	adjust_pvalue(method = "bonferroni") %>%  # 多重检验校正
	add_significance()           # 添加显著性标记
stat.test
  • group_by(dose): 在每个剂量水平下分别进行分析
  • t_test(len ~ supp): 执行t检验，比较VC和OJ两种补充类型
  • adjust_pvalue(): 使用Bonferroni法控制多重比较错误
  • 输出包含：统计量、自由度、p值、校正后p值、显著性标记

基础可视化

# Create a box plot
bxp <- ggboxplot(
	df, x = "supp", y = "len", fill = "#00AFBB", 
	facet.by = "dose"  # 按剂量分面显示
)
 
# Make facet and add p-values
stat.test <- stat.test %>% add_xy_position(x = "supp")
bxp + stat_pvalue_manual(stat.test)
  • facet.by = "dose": 创建三个子图，分别对应0.5、1、2 mg剂量
  • add_xy_position(): 自动计算p值标注的位置
  • 蓝色箱线图清晰展示数据分布和组间差异

高级可视化定制

bxp <- ggboxplot(
  df, x = "supp", y = "len", fill = "#00AFBB", 
  facet.by = "dose", scales = "free", add = "jitter"
  )
bxp +  
  stat_pvalue_manual(
    stat.test, bracket.nudge.y = -2, hide.ns = TRUE,
    label = "{p.adj}{p.adj.signif}"
    ) +
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))
  • scales = "free": 每个分面使用独立的y轴刻度
  • add = "jitter": 添加抖动点显示原始数据分布
  • bracket.nudge.y = -2: 调整括号的垂直位置
  • hide.ns = TRUE: 隐藏不显著的结果
  • label = "{p.adj}{p.adj.signif}": 同时显示具体p值和显著性星号
  • expansion(): 调整y轴扩展范围，为p值标签留出空间

条形图展示

# Create a bar plot with error bars (mean +/- sd)
bp <- ggbarplot(
	df, x = "supp", y = "len", add = c("mean_sd", "jitter"), 
	fill = "#00AFBB", facet.by = "dose"
)
 
# Add p-values onto the bar plots
stat.test <- stat.test %>% add_xy_position(fun = "max", x = "supp")
bp + stat_pvalue_manual(stat.test)

条形图特点：

  • add = c("mean_sd", "jitter"): 显示均值±标准差和原始数据点
  • fun = "max": p值标签位置基于数据最大值计算
  • 适合展示均值和变异程度

3. 双分组变量的复杂分析

双重分组统计

stat.test <- df %>%
	group_by(group, dose) %>%    # 按两个变量分组
	t_test(len ~ supp) %>%       # 在每个组合内比较supp
	adjust_pvalue(method = "bonferroni") %>%
	add_significance()
stat.test

分析层次：

  • 第一层分组: group (grp1, grp2)
  • 第二层分组: dose (0.5, 1, 2)
  • 在每个子组内: 比较supp (OJ vs VC)

双重分组可视化

stat.test <- stat.test %>% add_xy_position(x = "supp")
ggboxplot(
	df, x = "supp", y = "len", fill = "#E7B800",
	facet = c("group", "dose")  # 双重分面
) +
	stat_pvalue_manual(stat.test, hide.ns = TRUE)

双重分面说明：

  • facet = c("group", "dose"): 创建2×3的图形矩阵
  • 行: group (grp1, grp2)
  • 列: dose (0.5, 1, 2)
  • 每个子图内比较两种supp类型

双重分面条形图

stat.test <- stat.test %>% add_xy_position(x = "supp", fun = "mean_sd")
ggbarplot(
	df, x = "supp", y = "len", fill = "#E7B800",
	add = c("mean_sd", "jitter"), facet = c("group", "dose")
) +
	stat_pvalue_manual(stat.test, hide.ns = TRUE)
  • fun = "mean_sd": 基于均值+标准差计算标签位置
  • 金色条形图与之前的蓝色系形成对比

4. 多分组的多重比较

多组间比较

stat.test <- df %>%
	group_by(supp) %>%          # 按补充类型分组
	t_test(len ~ dose)          # 在每个supp内比较三个剂量
stat.test

比较类型：

  • 在OJ组内: 比较三个剂量水平
  • 在VC组内: 比较三个剂量水平
  • 使用方差分析后的事后检验可能更合适

多组比较可视化

# Box plots with p-values
stat.test <- stat.test %>% add_y_position()
ggboxplot(df, x = "dose", y = "len", fill = "#FC4E07", facet.by = "supp") +
	stat_pvalue_manual(stat.test, label = "p.adj.signif", tip.length = 0.01) +
	scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))
  • add_y_position(): 为多组比较自动计算y轴位置
  • tip.length = 0.01: 调整显著性线的尖端长度
  • 橙色系填充，视觉突出

多组条形图

ggbarplot(
	df, x = "dose", y = "len", fill = "#FC4E07",
	add = "mean_sd", facet.by = "supp"
) + 
	stat_pvalue_manual(stat.test,  label = "p.adj.signif", tip.length = 0.01) +
	scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))

5. 复杂实验设计的完整分析

三重分组统计

stat.test <- df %>%
	group_by(group, supp) %>%    # 双重分组
	t_test(len ~ dose) %>%       # 三组比较
	adjust_pvalue(method = "bonferroni") %>%
	add_significance()
stat.test

分析复杂度：

  • 2个group × 2个supp = 4个分析组
  • 每个组内进行3个剂量水平的比较

三重分组可视化

stat.test <- stat.test %>% add_xy_position(x = "dose")
ggboxplot(
	df, x = "dose", y = "len", fill = "#FC4E07",
	facet = c("group", "supp")   # 2×2分面矩阵
) +
	stat_pvalue_manual(stat.test, hide.ns = TRUE) +
	scale_y_continuous(expand = expansion(mult = c(0.05, 0.15)))

关系图谱

反向链接