浅析R语言多组定量资料非参数检验的多组比较及簇状柱形图显著性字母标记之分面与分组的图形艺术

R语言多组定量资料非参数检验的多组比较

非参数检验的应用

本流程是在刘永鑫老师提供的代码资料指导下完成

测试数据及代码链接：https://pan.baidu.com/s/1s0IX2eVe34jL2DViWCvZfQ
提取码：dunw

链接失效，后台回复"多组比较"，获得更新链接。

先简单说明一下非参数检验

参数检验是在已知总体分布的条件下(一般要求总体服从正态分布)对一些主要的参数(如 均值、百分数、方差、相关系数等)进行的检验，有时还要求某些总体参数满足一定的条件。如独立样本的t检验和方差分析不仅要求总体符合正态分布，还要求各总体方差齐性。

而当我们的数据不满足参数检验时我们可以选择非参数检验-秩和检验(Wilcoxon/Kruskal-Wallis H)

非参数检验方法简便，不依赖于总体分布的具体形式，因而 适用性强，但灵敏度和精确度不如参数检验。一般而言，非参数检验适用于以下3种情况

(1)顺序类型的数据资料，这类数据的分布形态一般是未知的。

(2)虽然是连续数据，但总体分布形态未知或者非正态，这和卡方检验一样，称为自由分布检验。

(3)总体分布虽然正态，数据也是连续类型，但样本容量极小，如10以下(虽然T检验被
这叫小样本统计法，但样本量太小，代表性毕竟很差，最好不要用比较严格的参数检验法)。

非参数检验的缺点与不足

由于这些特点，再加上非参数检验法的一般原理和计算方法比较简单，所以常被用于数据统计的一些前期研究中进行正式研究。当然，由于许多非参数检验不涉及参数计算，数据中的信息没有得到充分利用，因此其统计检验能力远不如参数检验。

关于使用非参数检验的纠结

在很多人看来，参数检验似乎是统计方法的"主流"， 而非参数检验往往被当成"非主流"。大家似乎更喜欢用t检验、方差分析这样的参数检验。 即使在数据不满足正态分布的条件下，仍然有使用参数检验的执念；总觉得非参数检验不够正式，或者就是不喜欢用。

而现实中的很多数据都 不服从正态分布，这时候 非参数检验的统计学效能要高于参数检验。并且， 非参数检验更加稳健。参数统计方法是建立在严格假设条件基础上，一旦假设条件不符合，其推断的正确性就不存在了。 非参数检验带有最弱的假设，对模型的限制很少，因而天然的具有稳健性。

多组定量资料非参数检验及多组比较的实现

kruskal.test()函数可以确定总体差异是否有统计学意义，但是不知道具体那些组之间存在差异。可以使用wilcox.test()函数每次进行两组数据比较。 一种更为有效的方法是在控制犯第一类错误的概率前提下,执行可以同步进行的多组比较,这样可以直接完成所有组之间的成对比较。 nparcomp包提供了所需要的非参数多组比较程序。此包中的 nparcomp()函数接受的输入为一个两列的数据框,其中一列为因变量,另一列为分组变量。 ==本推文的第三种方法采用此方法==。

本文分别采用三个统计R包的三种不同函数方法进行非参数多组比较：

• "pgirmess"统计包的kruskalmc()函数
• "PMCMRplus"统计包的posthoc.kruskal.nemenyi.test()函数
• "nparcomp"统计包的nparcomp()函数

自己挑选喜欢的用

library(tidyverse)#数据整理与数据转换用了一些更好用更易懂的函数
library(ggprism)
dat %
  dplyr::select(SampleID,Group) %>%
  inner_join(dat,by='SampleID')#按照SampleID内连接文件
head(dat, n = 3)
dat

函数和参数简介

函数参数设置：

• data就是上面整好的数据,
• group是你的分组信息列，比如α多样性的种类(或不同的基因),
• compare是每个α多样性要比较的不同处理(或每个gene要比较的不同处理),
• value 值就是要比较的α多样性/gene拷贝数的数值。

整体思想如下(例如本数据)：
首先给输入数据dat，根据alpha列分成不同的小子集，每个小子集比较不同Group下v值的差异情况，最后汇总输出。

"pgirmess"统计包的kruskalmc()函数

# 1 -----------------------------------------------------------------------
kruskalmc_compare1

"PMCMRplus"统计包的posthoc.kruskal.nemenyi.test()函数

# 2 -----------------------------------------------------------------------
PMCMR_compare1 %
      mutate(compare=rownames(n)) %>%
      gather(group,p,-compare,na.rm = TRUE) %>%
      unite(compare,group,col="G",sep="-")
    dif

"nparcomp"统计包的nparcomp()函数

# 3 -----------------------------------------------------------------------
nparcomp_compare1 %
        Replace(from = ' , ',to='-') %>%
        Replace(from=c('p\\(','\\)'),to='')#正则表达式替换
    temp_name %
                 separate(col = 'tn',sep = '-',into = c('n1','n2'))
    difname %
                Replace(from = ' - ',to='-')#正则表达式替换
    names(dif)

alpha多样性在不同处理下的差别

运行，这里拿alpha多样性测试，看不同alpha多样性在不同处理下的差别。

#df1

在此可以查看各个α多样性下不同处理间的多组比较字母标注结果，这也是本脚本的亮点之一

数据量很大的情况下，可以直接查看差异情况，不用一个个的做出图再点开看，很是方便。

p1 = ggplot(dat)+geom_boxplot(aes(x=Group,y=v,fill=Group))+
  geom_text(data=df1,aes(x=Group,y=mean+1.3*std,label=Letters))+
  facet_wrap(~alpha,scales = "free_y")+ labs(x='Group',y='AlphaDiv')+
  ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p1

本图一张即可包含所有数据情况，方便查看

#df2 b>c)

p2 = ggplot(dat)+geom_boxplot(aes(x=Group,y=v,fill=Group))+
  geom_text(data=df2,aes(x=Group,y=mean+1.3*std,label=Letters))+
  facet_wrap(~alpha,scales = "free_y")+ labs(x='Group',y='AlphaDiv')+
  ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p2

#df3 b>c)

p3 = ggplot(dat)+geom_boxplot(aes(x=Group,y=v,fill=Group))+
  geom_text(data=df3,aes(x=Group,y=mean+1.3*std,label=Letters))+
  facet_wrap(~alpha,scales = "free_y")+ labs(x='Group',y='AlphaDiv')+
  ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p3

簇状柱形图显著性字母标记之分面与分组的图形艺术

只是展示方法(多样性数据可能不适合做这个图，但是纵轴标度一致的数据且有需求要做，适合做这个图)。

颜色设置

library(RColorBrewer)
library("scales")
#自己选颜色
display.brewer.all()

#选取颜色
display.brewer.pal(9,"Set1")

#可以直接选取
yanse=brewer.pal(9,"Set1")
#也可以显示出色号后
show_col(brewer.pal(9,"Set1"))

#再挑选喜欢的色号
yanse=c("#E41A1C","#377EB8","#4DAF4A",
     "#984EA3","#FF7F00","#FFFF33",
     "#A65628","#F781BF","#999999")

分面的图形艺术

数据用的df2

df1$alpha

分组的图形艺术

数据用的df2

p5 = ggplot(df2,aes(x =alpha,y = mean,fill = Group))+
  geom_bar(stat = "identity", width=0.9, position = position_dodge(0.9)) +
  geom_text(aes(label = Letters),size=7,vjust=-1.2,position = position_dodge(0.9))+
  geom_errorbar(aes(ymin = mean - std/3^0.5,ymax = mean+std/3^0.5),
                position = position_dodge(0.9), width = .1,color="red",lwd=1)+
  labs(x='Group',y='AlphaDiv')+
  scale_y_continuous(expand = c(0,0),limits = c(0,H))+
  scale_fill_manual(values = yanse, guide = guide_legend(title = NULL))+
  ggprism::theme_prism()+theme(axis.text.x = element_text(angle = 45))
p5

Output figure width and height
Letter纸图片尺寸为单栏89 mm，双栏183 mm，页面最宽为247 mm
推荐比例16：10，即半版89 mm x 56 mm; 183 mm x 114 mm

# ggsave("./alpha1.pdf", p1, width = 350, height = 200, units = "mm")
# ggsave("./alpha2.pdf", p2, width = 350, height = 200, units = "mm")
# ggsave("./alpha3.pdf", p3, width = 350, height = 200, units = "mm")
# ggsave("./alpha4.pdf", p4, width = 600, height = 200, units = "mm")
# ggsave("./alpha5.pdf", p5, width = 600, height = 200, units = "mm")

猜你想学

R统计-正态性分布检验[Translation]

R统计-方差齐性检验[Translation]

参考资料

《R语言统计分析与应用》

浅析R语言单因素方差分析中的多重比较

在R语言上做Kruskal–Wallis秩和检验

Wilcoxon, Friedman...扒一扒非参数检验名称中的牛人

作者：flyfly、旭日阳光

责编：马腾飞 南京农业大学

审核：刘永鑫 中科院遗传发育所

猜你喜欢

10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature Cell专刊 肠道指挥大脑

系列教程：微生物组入门 Biostar 微生物组 宏基因组

专业技能：学术图表 高分文章 生信宝典 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索 Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun

在线工具：16S预测培养基 生信绘图

科研经验：云笔记 云协作 公众号

编程模板: Shell R Perl

生物科普: 肠道细菌 人体上的生命 生命大跃进 细胞暗战 人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了"宏基因组"专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注"姓名-单位-研究方向-职称/年级"。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注"宏基因组"

点击阅读原文，跳转最新文章目录阅读

Original: https://blog.csdn.net/woodcorpse/article/details/120232231
Author: 刘永鑫Adam
Title: 浅析R语言非参数检验的多组比较及分面与分组的图形艺术

相关阅读

Title: 图像分类篇——AlexNet详解

一、概述

AlexNet是由2012年ImageNet竞赛参赛者Hinton和他的学生Alex Krizhevsky设计的。AlexNet在当年赢得了ImageNet图像分类竞赛的冠军，使得CNN成为图像分类问题的核心算法模型，同时引发了神经网络的应用热潮。
1. AlexNet的创新
作为具有历史意义的网络结构，AlexNet包含以下方面的创新。
(1)非线性激活函数 ReLU
在AlexNet出现之前，sigmoid是最为常用的非线性激活函数。sigmoid函数能够把输入的连续实值压缩到0和1之间。但是，它的缺点也非常明显：当输入值非常大或者非常小的时候会出现饱和现象，即这些神经元的梯度接近0，因此存在梯度消失问题。为了解决这个问题，AlexNet使用ReLU作为激活函数。
ReLU函数的表达式为F(x)=max(0，z)。若输入小于0，那么输出为0；若输入大于0，那么输出等于输入。由于导数始终是1，会使得计算量有所减少，且AlexNet的作者在实验中证明了，ReLU 函数的收敛速度要比sigmoid 函数和 tanh函数快。
(2)局部响应归一化
局部响应归一化(local response normalization，LRN)的思想来源于生物学中的"侧抑制"，是指被激活的神经元抑制相邻的神经元。采用LRN的目的是，将数据分布调整到合理的范围内，便于计算处理，从而提高泛化能力。虽然ReLu函数对较大的值也有很好的处理效果，但AlexNet的作者仍然采用了LRN的方式。式3.7是Hinton在有关AlexNet的论文中给出的局部响应归一化公式。

a表示第i个卷积核在(x，y)处经卷积、池化、ReLU函数计算后的输出，相当于该卷积核提取的局部特征。N表示这一层的卷积核总数；n表示在同一位置的临近卷积核的个数，是预先设定的；k、alpha、beta均为超参数。假设N=20，超参数和n按照论文中的设定，分别为k=2、alpha=10-4、beta=0.75、n=5。第5个卷积核在(x，y)处提取了特征，的作用就是以第5个卷积核为中心，选取前后各5/2=2个(取整，因此n所指个数包含中心卷积核)卷积核，因此有max=3和min=7，卷积核个数就是3、4、5、6、7。
通过分析公式，我们对局部响应归一化的理解会进一步加深。式3.7的存在，使得每个局部特征最后都会被缩小(只是缩小的比例不同)，相当于进行了范围控制。一旦在某个卷积核周围提取的特征比它自己提取的特征的值大，那么该卷积核提取的特征就会被缩小。相反，如果在该卷积核周围提取的特征比它自己提取的特征的值小，那么该卷积核提取的特征被缩小的比例就会变小，最终的值与周围的卷积核提取的特征值相比就显得比较大了。
(3) dropout
        dropout通过设置好的概率随机将某个隐藏层神经元的输出设置为0，因此这个神经元将不参与前向传播和反向传播，在下一次迭代中会根据概率重新将某个神经元的输出置0。这样一来，在每次迭代中都能够尝试不同的网络结构，通过组合多个模型的方式有效地减少过拟合。
(4)多GPU训练
        单个GPU的内存限制了网络的训练规模。采用多GPU协同训练，可以大大提高AlexNet的训练速度。

二、AlexNet的结构

下面介绍AlexNet的结构。如图所示，由于采用双GPU协同训练，该网络结构图分为上下两部分，且两个GPU只在特定的层内通信。该模型一共分为8层，包括5个卷积层和3个全连接层，每个卷积层都包含激活函数ReLU、池化和LRN处理。

AlexNet架构 层类型特征图大小内核步幅填充激活OUT全连接-1000---ReLUF10全连接-4096---ReLUF9全连接-4096---ReLUs8最大池化25666332valid-C7卷积2561313551sameReLUC6卷积3841313551sameReLUC5卷积3841313551sameReLUS4最大池化2561313332valid-C3卷积2562727551sameReLUS2最大池化962727332valid-C1卷积96555511114validReLUIn输入3(RGB)227*227----

(1)Layer1卷积层
该层接收227×227×3(R、G、B三个通道)的输入图像，使用96个大小为11×11×3的卷积核进行特征提取，步长为1，扩充值为0，通过式3.1可以得到输出特征图的尺寸。由于使用了双GPU，图3.13中的每个GPU处理的卷积核的数量均为48个。卷积后，得到大小为55×55×96的特征图。
对得到的特征图，使用ReLU函数将它的值限定在合适的范围内，然后使用3×3的滤波器进行步长为2的池化操作，得到27×27×96的特征图。最后，进行归一化处理，规模不变。
(2)Layer 2卷积层
Layer 2与Layer 1相似，它接收Layer 1输出的27×27×96的特征图，采用256个5×5×96的卷积核，步长为1，扩充值为2，边缘用0填充，生成27×27×256的特征图。紧随其后的是与Layer 1相同的ReLU函数和池化操作，最后进行归一化处理。
(3) Layer3卷积层与Layer 4卷积层
Layer3和Layer4也采用了卷积层，但只进行卷积和ReLU。不进行池化和归一化操作。Layer3的每个GPU都有192个卷积核，每个卷积核的尺寸是3×3×256.步长为1，扩充值为1，边缘用0填充。最终，每个GPU都生成的个13×13的特征图。
Layer 4与Layer3的区别在于卷积核的尺寸。Layer4不像Layer3那样接收的一层所有GPU的输入，而只接收所在GPU的输出。因此，Layer 4的卷积核的尺寸为3×3×192，每个GPU都有192个卷积核。
与Layer 3相同的是，Layer 4仍然进行扩充值为1的0填充，且步长为1。最终，Layer 4的每个GPU都生成192个13×13的特征图。卷积后，这两个层都会由ReLU函数进行处理。
(4)Layer 5卷积层
Layer 5会依次进行卷积、ReLU和池化操作，但不进行归一化操作。该层中的每个GPU都接收本GPU中 Layer 4的输出，每个GPU使用128个3×3×192的卷积核，步长为1，使用扩充值为1的0填充，各自生成128个13×13的特征图后进行池化操作。池化尺寸为3×3，步长为2，最终生成6×6×128的特征图(两个GPU，共256个)。
(5)Layer 6全连接层
从该层开始，后面的网络层均为全连接层。Layer 6仍然按GPU进行卷积，每个GPU使用2048个6×6×256的卷积核，这意味着该层中的每个GPU都会接收前一层中两个GPU 的输出。卷积后，每个GPU都会生成2048个1×1的特征图。最后，进行ReLU和dropout操作，两个GPU共输出4096个值。
(6)Layer 7全连接层
Layer7与Layer6相似，与Layer6进行全连接，在进行ReLU和dropout操作后，共输出4096个值。
(7)Layer 8全连接层
Layer 8只进行全连接操作，且该层拥有1000个神经元，最终输出1000个float型的值(该值即为预测结果)。

三、构建AlexNet

构建LRN层

LRN层通常用lambda层进行包装，参数一般设置如下：

model.add(keras.layers.lambdas(lambda x: tf.nn.lrn(x, depth_radius=2, alpha=0.00002, beta=0.75, bias=1)))

也可以用自定义层实现：

class LRN(keras.layers.Layer):
    def __init__(self, depth_radius=2, bias=1, alpha=0.00002, beta=0.75, **kwargs):
        super().__init__(**kwargs)
        self.depth_radius = depth_radius
        self.bias = bias
        self.alpha = alpha
        self.beta = beta

    def call(self, input):
        return tf.nn.lrn(input, self.depth_radius, self.bias, self.alpha, self.beta)

    def get_config(self):
        base_config = super().get_config()
        return {**base_config, 'depth_radius': self.depth_radius,
                'bias': self.bias, 'alpha': self.alpha, 'beta': self.beta}

构建模型：

model = keras.Sequential()
# 第一层
model.add(keras.layers.Conv2D(filters=96, kernel_size=(11,11), strides=(4,4), input_shape=(227,227,3), padding="VALID", activation="relu"))
model.add(LRN())
model.add(keras.layers.MaxPooling2D(pool_size=(3,3), strides=(2,2), padding="same"))
# 第二层
model.add(keras.layers.Conv2D(filters=256, kernel_size=(5,5), strides=(1,1), padding="SAME", activation="relu"))
model.add(LRN())
model.add(keras.layers.MaxPooling2D(pool_size=(3,3), strides=(2,2), padding="same"))
# 第三层
model.add(keras.layers.Conv2D(filters=384, kernel_size=(5,5), strides=(1,1), padding="SAME", activation="relu"))
# 第四层
model.add(keras.layers.Conv2D(filters=384, kernel_size=(5,5), strides=(1,1), padding="SAME", activation="relu"))
# 第五层
model.add(keras.layers.Conv2D(filters=256, kernel_size=(5,5), strides=(1,1), padding="SAME", activation="relu"))
model.add(keras.layers.MaxPooling2D(pool_size=(3,3), strides=(2,2), padding="same"))
# 第六层
model.add(keras.layers.Flatten())
model.add(keras.layers.Dense(4096, activation="relu"))
model.add(keras.layers.Dropout(0.5))
# 第七层
model.add(keras.layers.Dense(4096, activation="relu"))
model.add(keras.layers.Dropout(0.5))
# 第八层
model.add(keras.layers.Dense(10, activation="softmax"))

注：LRN层也可用BatchNormalization替代，性能也还OK

Original: https://blog.csdn.net/qq_52053775/article/details/124377926
Author: 樱花的浪漫
Title: 图像分类篇——AlexNet详解