前言
上一期我们说：60-R可视化-8-用ggsignif做统计分析绘图
对于已有的原始数据进行绘图非常的方便。
可是，如果我们拿到手的就是处理后的统计结果呢？
这时候需要我们自己计算一下了。
正常操作
手动计算t值与p值
比如通过REdaS::freqCI 方法获得的结果：
代码语言：javascript复制> 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等。
我的数据如下：
代码语言：javascript复制> 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 检验的计算公式：
完整代码如下：
代码语言：javascript复制# 解决粉丝儿的一个小问题
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 竟然不提供外部统计结果的调整位置函数。手动调整太痛苦了：
代码语言：javascript复制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))
另一个花招
既然数据是符合正态的，那我们为每个检验的数据，生成若干个仿真数据，用仿真数据来进行绘图检验，从理论上来看，也是可以的。
至于这个若干个数据数值设定为多少，还需要具体考虑这个统计结果来自何种分布的数据，具体问题具体对待。
至于本例中，freqCI 其实就是从正态抽取了频数个数的数据，那我们将数值设置为相同的频率个数N即可，那么自由度也就是N-1。
先挖个坑~
我的思考
ggsignif 虽然没有给出它实现绘图统计显著注释棒自动调整函数的接口，但实际上我们或许可以通过它的源代码，来实现自己计算的统计结果绘图的自动调整。
再挖个坑~
拓展阅读
关于在R中手动计算p值与 t值：用R语言解读统计检验-T检验 | 粉丝日志 (fens.me)粉丝日志 (fens.me)[1]
参考资料
[1]粉丝日志 (fens.me): http://blog.fens.me/r-test-p
