在仙桃学术上的生信工具里面，有一个折线图的绘图工具，可以很快速便捷的得出结论并可视化结果，当然不是说这个功能有多强大，而是统计学方法非常专业。
比如用它自带的数据https://bioinfomatics.xiantao.love/biotools/data/demo/free/linePlot/%E6%8A%98%E7%BA%BF%E5%9B%BE.xlsx
通过无脑式的鼠标点击，可得到下面一系列的结果。

折线图

比如统计描述

又比如异常值分析、正态性检验(Shapiro-Wilk normality test)、方差齐性检验(Levene's test)、协方差同质假设(Homogeneity of covariances assumption/Box’s M-test)、球形假设(Mauchly's Test for Sphericity)、两因素单向重复测量数据方差分析（Two-way mixed ANOVA）、单独效应分析(Simple effect)、多重成对比较(Pairwise Comparisons of Estimated Marginal Means)等等专业术语，并且还给了详细的统计学结果。

image.png

image.png

image.png

这种统计学结果让大家很欣慰，那么这些结果都是如何计算出来的呢？？？

我们可以点开demo.R这个结果，读一读代码
https://bioinfomatics.xiantao.love/biotools/code/open/lineplot.R

现在我们复现一下：

加载需要的包并导入数据

library(ggplot2)
library(reshape2)
library(car)
library(rstatix)
## 导入数据，比如我们把下载的折线图数据保存在桌面
library(readxl)
data <- read_excel("~/Desktop/折线图.xlsx") ## mac系统代码
data # 显示结果

新增一列id，id即为数字，用于后续分析，必不可少

data$id <- 1:nrow(data)

将短数据转换为长数据

data2 <- gather(data, key = "x", value = "value", -trt, -id)

对各组进行统计学描述

data3 <- data2 %>% 
  group_by(trt, x) %>% 
  get_summary_stats(value)

批量运行正态性检验（Shapiro-Wilk normality test）

for(i in unique(data[,1])){
  data1 <- data[data[,1] == i,]
  print(lapply(data1[,-1], function(x) shapiro.test(x)))
}

$Time1

Shapiro-Wilk normality test

data: x

W = 0.91913, p-value = 0.2788

$Time2

Shapiro-Wilk normality test

data: x

W = 0.88011, p-value = 0.08792

$Time3

Shapiro-Wilk normality test

data: x

W = 0.73897, p-value = 0.002055

$id

Shapiro-Wilk normality test

data: x

W = 0.95933, p-value = 0.7742

批量运行方差齐性检验（Levene's Test）

for(i in unique(data2$x)){
  data1 <- data2[data2$x == i,]
  print(leveneTest(value~trt, data = data1))
}

Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 1.5617 0.2245
22
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 2.5864 0.122
22
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 3.8375 0.06291 .
22

---

Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1

两因素重复测量方差分析

anova_test(data = data2, dv = value, wid = id, within = x, between = trt) 

ANOVA Table (type II tests)

$ANOVA
Effect DFn DFd F p p<.05 ges
1 trt 1 22 510.657 1.02e-16 * 0.918
2 x 2 44 391.826 9.19e-29 * 0.902
3 trt:x 2 44 407.706 4.01e-29 * 0.905

$Mauchly's Test for Sphericity
Effect W p p<.05
1 x 0.966 0.699
2 trt:x 0.966 0.699

$Sphericity Corrections
Effect GGe DF[GG] p[GG] p[GG]<.05 HFe DF[HF] p[HF] p[HF]<.05
1 x 0.968 1.94, 42.57 6.65e-28 * 1.059 2.12, 46.61 9.19e-29 *
2 trt:x 0.968 1.94, 42.57 2.98e-28 * 1.059 2.12, 46.61 4.01e-29 *

可视化结果

ggplot(data = data3, aes(x = x, y = mean, color = trt, group = trt)) +
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), color = "black", width = 0.2) +
  geom_line() +
  geom_point() +
  theme_bw()

image.png

当然这个图跟原图不一样，并且没有两两比较，我们可以使用ggpubr和rstatix进行二次绘图。

关于rstatix进行统计学分析，其实可以Repeated Measures ANOVA in R: The Ultimate Guide - Datanovia
得到答案，具体的分析，我们暂不说了，只说如何进行统计学分析，并且自动两两比较

res.aov <- anova_test(data = data2, dv = value, wid = id, within = x, between = trt) 
# 组间两两比较
pwc <- data2 %>%
    group_by(x) %>%
    pairwise_t_test(
        value ~ trt, paired = TRUE,
        p.adjust.method = "bonferroni" ## 校正方法，包括 "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none".等
    )
pwc

可以看到已经自动进行了两两比较，我们使用rstatix包的add_xy_position()函数可以添加两两比较列表和x轴y轴位置。

 pwc <- pwc %>% add_xy_position(x = "x")# 按时间分列
pwc

可视化绘图

library(ggpubr)
ggline(data2,x='x',y='value',
color = 'trt', #按组配色
add = 'mean_sd', #添加均数标准差，也可以设置均数标准误，CI等。
palette = "aaas" # aaas杂志配色
)+stat_pvalue_manual(pwc, tip.length = 0, hide.ns = TRUE) + # 添加两两比较，隐藏无意义
    labs(
        subtitle = get_test_label(res.aov, detailed = TRUE), # 添加整体差异
        caption = get_pwc_label(pwc) # 右下角条件两两比较方法。
    )

image.png

当然，我们还可以继续修改图片，比如全部显示结果，取消显著性标记的下划线，或者显示具体的p值等

ggline(data2,x='x',y='value',
       color = 'trt', #按组配色
       add = 'mean_sd', #添加均数标准差，也可以设置均数标准误，CI等。
       palette = "aaas",# aaas杂志配色
ggtheme = theme_bw(), #背景
xlab = F,ylab = "Score", #xy轴名称
legend.title="Treatment" ,legend="top", #标签名称和位置
title = "Comparison of two groups at different times" #标题
)+stat_pvalue_manual(pwc, bracket.size = 0) + # 添加两两比较，显示所有结果
    labs(
        subtitle = get_test_label(res.aov, detailed = TRUE), # 添加整体差异
        caption = get_pwc_label(pwc) # 右下角显示事后两两比较方法。
    )

image.png

也可以显示数值

ggline(data2,x='x',y='value',
       color = 'trt', #按组配色
       add = 'mean_sd', #添加均数标准差，也可以设置均数标准误，CI等。
       palette = "aaas",# aaas杂志配色#背景
       xlab = F,ylab = "Score", #xy轴名称
       legend.title="Treatment" ,legend="right",
       title = "Comparison of two groups at different times" #标题
)+stat_pvalue_manual(pwc, bracket.size = 0,label = "p.adj") + # 添加两两比较，显示所有结果
    labs(
        subtitle = get_test_label(res.aov, detailed = TRUE), # 添加整体差异
        caption = get_pwc_label(pwc) # 右下角显示事后两两比较方法。
    )

image.png