61-R可视化-9-对已有统计结果数据做统计分析绘图

北野茶缸子

发布于 2021-12-17 10:57:11

3350

发布于 2021-12-17 10:57:11

文章被收录于专栏：北野茶缸子的专栏

前言

上一期我们说：60-R可视化-8-用ggsignif做统计分析绘图

对于已有的原始数据进行绘图非常的方便。

可是，如果我们拿到手的就是处理后的统计结果呢？

这时候需要我们自己计算一下了。

正常操作

手动计算t值与p值

比如通过REdaS::freqCI 方法获得的结果：

> tmp <- sample(c("a","b","c"), replace = T, size = 100)
> table(tmp)
tmp
 a  b  c 
37 35 28 
> tmp2 <- REdaS::freqCI(tmp, level = .95)
Warning: "x" was supplied as a "character" vector. A "factor"-type variable was assumed.

> tmp2
      2.5% Estimate    97.5%
a    27.54       37    46.46
b    25.65       35    44.35
c    19.20       28    36.80

比如这个结果，返回的就是95% 的置信区间。这里我们暂时不去看freqCI 函数算法本身是否正确。

通过这个结果我们可以得到什么呢？得到均值，通过95% 置信区间，是可以算出来标准差的。那么也就可以进行两组正态数据的假设检验，看A数据与B 数据均值是否相等。

方法选择就是常规操作。因为这里我们的数据是正态分布的，两个样本，就使用t 检验。如果是多组的话，那就去用anova等。

我的数据如下：

> head(my_tmp)
                   lower   Estimate     upper   Both          cell
Ccr6_ILC3      0.5294785  0.7902461  1.051014 female     Ccr6_ILC3
Foxp3_Treg    15.7633224 16.8661097 17.968897 female    Foxp3_Treg
ICL2           4.1155780  4.7414766  5.367375 female          ICL2
Il17_ILC3s    18.5173287 19.6884172 20.859506 female    Il17_ILC3s
LTi_like_ILC3 13.0862402 14.1115376 15.136835 female LTi_like_ILC3
NK             4.7095683  5.3736735  6.037779 female            NK

我实际需要进行假设检验的分组为，每个cell 组别下的不同Both 列之间的数据进行比较：

这里我们首先看看两独立样本t 检验的计算公式：

完整代码如下：

# 解决粉丝儿的一个小问题
load("./to_pp_from_zzc.Rdata")
tab.1$Both <- as.character(tab.1$Both)
my_tmp <- tab.1[,c(-6)]
# my_tmp <- as.tibble(my_tmp)
# my_tmp$df <- my_tmp[,1:3]
# 这里还需要一下每个数据的数目
# 因为我们是通过频率计算的统计结果，这里也就是各个分组的频数
# 粉丝儿并没有给我结果
# 这里我统一设为50
my_tmp$freq <- 50
my_list <- split(my_tmp, list(my_tmp$cell))

my_T.Test <- function(a1, a2, a3, a4) {
  # a1 is a vector containing mean lower upper
  # so does a2, other group
  mean1 <- a1[2]
  mean2 <- a2[2]
  lower1 <- a1[1]; lower2 <- a2[1]
  upper1 <- a1[3]; upper2 <- a2[3]
  thita1 <- (upper1 - lower1)/ 2# betweeness of 95% lower and upper is equal to 2*sd 
  thita2 <- (upper2 - lower2)/ 2
  t <- as.numeric((mean1 - mean2)/(sqrt((thita1**2/a3) + (thita2**2/a4)))) # 计算t值
  df <- a3+a4-2 # 计算自由度
  # p <- 2*pt(-abs(t_stat),df=n-1)
  p <- pt(-abs(t), df = df)  # pt 函数计算p值
}


my_re1 <- lapply(my_list, function(x){
  group_name <- x[,"cell"]
  a1 <- as.numeric(x[1,1:3])
  a2 <- as.numeric(x[2,1:3])
  a3 <- x[1,"freq"]
  a4 <- x[2,"freq"]
  p <- round(my_T.Test(a1,a2,a3,a4),5)
})

my_re2 <- do.call("rbind", my_re1)、

> head(my_re2)
                 [,1]
Ccr6_ILC3     0.05060
Foxp3_Treg    0.00000
ICL2          0.00024
Il17_ILC3s    0.00006
LTi_like_ILC3 0.00000
NK            0.00000

现在获得了p值，我们就可以添加到图表上了。

手动添加：痛苦源泉

上一讲我们提到ggsignif 是可以支持手动添加的。

所以，我们现在就得把这个annotations 给手动加上去了！

这里非常鸡肋啊，ggsignif 竟然不提供外部统计结果的调整位置函数。手动调整太痛苦了：

my_re1 <- lapply(my_list, function(x){
  group_name <- x[,"cell"]
  a1 <- as.numeric(x[1,1:3])
  a2 <- as.numeric(x[2,1:3])
  a3 <- x[1,"freq"]
  a4 <- x[2,"freq"]
  p <- round(my_T.Test(a1,a2,a3,a4),5)
})

my_re2 <- do.call("rbind", my_re1)
head(my_re2)

ggplot(tab.1, aes(cell, Estimate)) + geom_boxplot(aes(color = Both), position = "dodge") +
  geom_errorbar(aes(ymin=lower, ymax=upper, group = Both), width=.2,position=position_dodge(.7)) +
  geom_signif(annotations = my_re2, y_position = rep(23, 11), xmin = seq(0.7,10.7, 1), 
              xmax = seq(1.1,11.1, 1))