5.1假设检验和功效分析
①假设检验：统计学的主要目的是根据样本信息对未知的总体参数进行推断：从观察到的样本数据中得出关于总体的结论。统计推断包括总体参数的估计和关于总体参数的假设检验。假设检验是一种用来检验统计假设的统计程序。这是制定统计程序以决定接受或拒绝假设的统计推理的主要领域，而幂和样本量的计算与假设检验有关。统计假设是关于人口的断言或理论，特别是关于人口分布的参数，如位置(均值)、尺度(方差/离散度)。本章的目的是在一般情况下，特别是在微生物群研究中，介绍假设检验和相关的功率和样本量计算。
统计测试的要素：统计检验有四个要素：(I)零假设，H0；(Ii)替代假设，H1；(Iii)检验统计，和(Iv)拒绝区域。每个假设检验程序由两个假设组成。第一个假设是零假设(用H0表示)，这是一种关于一个或多个总体参数的特定值的理论。理论通常表述为H0：参数=特定值，例如H0：l？0：25。我们想要根据样本中包含的信息来检验的假设，通常是在希望它被拒绝的情况下制定的。第二个假设被称为替代假设(用Ha表示)，这是一种通常以与零假设相反的形式提出的理论。通常，零假设声称没有差异或没有变化，而另一种假设则声明有一些差异或变化已经发生。因此，零假设总是声明总体参数与索赔值完全相等，而另一种假设允许该参数有几个不同的值。例如，在一项临床试验中，如果研究人员想知道与接受安慰剂20 mm的患者相比，服用新药的患者是否降低了平均标准化平均得分的0.20分以下。零假设将是H0：l？0：20，而替代假设将是以下之一：HA：l[0：20，其中l表示总体的平均标准化平均分数。检验统计量是用于决定是否拒绝零假设的样本统计量。它是样本观测值和一些已知常数的函数，统计决策将基于这些常数。拒绝区域指定将拒绝零假设的测试统计量的数值。具体地说，它是一个概率区域，它告诉我们我们的理论或假设可能是真的，也可能不是真的。那么问题是：如何选择拒绝区域？简而言之，排斥区域与α水平相关。您需要选择您愿意接受的Alpha级别。例如，如果您想要95%确信您的结果是重要的，您将选择5%的alpha级别(100%-95%)。“5%水平”是拒绝区域。
检验统计假设的步骤：从样本统计到总体参数估计，统计假设检验起着基础性的作用。您可以使用两种方法进行假设检验：使用临界值和使用p值。对于临界值，推理步骤如下：(1)由于初始研究假设未知，说明相关的无效假设和备选假设；(2)对测试样本进行统计假设，如统计独立性或观测值的分布；(3)选择合适的检验，并指定相应的检验统计量T；(4)从假设出发，推导出零假设下检验统计量的分布；(5)选择通常为小值(例如0.01、0.05)的显著性水平(临界区)a，其被称为检验的显著性水平或低于其将拒绝零假设的概率阈值；(6)由观测值计算检验统计量T的观测值；(7)决定要么拒绝零假设，转而支持备选方案，要么接受它。您还可以决定拒绝或接受具有p值的零假设。如果要使用p值方法，则需要计算p值。这是在零假设下，获得等于或比实际观察到的结果更极端的结果的概率。如果p值落在拒绝区域，则意味着您具有统计上显著的结果；因此，您可以拒绝零假设，并得出结论，认为替代假设为真。如果p值落在拒绝区域之外，则意味着您的结果不能提供足够的证据来反对无效假设。
微生物组数据中的假设检验：对于微生物群落比较，一般假设测试可以写成：(1) H0: π1=π0 versus HA: π1不等于π0
(2) H0: π1=...πj=...πJ versus HA: πi不等于πj。上述两种假设检验分别类似于经典统计学中的单样本t检验和两样本t检验或方差分析，它们构成了微生物群研究中用于比较的基本统计假设检验框架。A.分类群与先前指定的微生物群落的平均比例 B.两个样本集分类群的平均比例 C.来自两个以上类群的平均分类群比例 D.分类群的平均比例和比例都是跨组的。
微生物组研究中的功率和样本量计算基于假设检验框架。微生物组数据既有其独特性，又有其他研究领域的共性。因此，假设检验和幂以及样本量计算具有相似的设置。对于微生物组数据的共同特征，根据这些数据值的分布方式和要比较的组数，可以对微生物组假设使用标准t检验、方差分析(ANOVA)或相应的非参数检验。对于微生物组数据的独特特征，研究人员试图开发适当的统计分析工具，包括功率和大小计算，以更好地拟合数据。例如，考虑到微生物组数据结构，开发了多变量分析工具来考虑分类群之间的相互作用或相互关系。其中，HMP软件包采用基于分类群计数的Dirichlet-多项式模型的参数方法。
②功效分析与样本量计算
功效分析与样本量计算的重要性：功率分析最常见的目的是确定合理检测给定大小的影响所需的最小对象。此外，如果有可用的效果大小和受试者数量，则使用功率分析来确定功率。功率分析还可以用来比较不同的统计检验过程，例如，在同一假设的参数检验和非参数检验之间。例如，当你为NIH拨款提案计划一项研究或设计一项实验时，你会问一个常见的问题：“我们需要多少科目？”这个问题很重要，因为样本大小应该足够大，如果存在科学兴趣的效应大小，它就有很好的机会被检测出来。以小样本进行研究，在时间和金钱上都是浪费资源，因为结果往往是没有定论的。但是，也不建议使用大样本量。首先，它不会有什么用处，因为在科学上重要性不大的效应大小在统计上也可能是有意义的。其次，出于经济原因，确定合适的样本量是很重要的。如果能从较小的样本中准确地找到答案，那就是浪费有限的可用资源。第三，在人体研究或随机对照试验中，招募比所需更多的受试者可能是不道德的。参与治疗的患者不应被误用，特别是安慰剂组患者。
功效分析：功效是当零假设(H0)实际上为假或替代假设(HA)为真时，测试正确拒绝该零假设(H0)的概率。它可以等效地被认为是当备选假设(HA)为真时接受该假设的概率；也就是说，能力是测试检测效果的能力，如果该效果确实存在的话。因此，要理解权力的概念，我们必须理解假设检验的原理。用于统计假设检验的四种可能的决策结果及其相关概率之间的关系如表5.1所示。
如果零假设H0是否真的为真，则组间无显著差异或变量间无关系，不应予以否定。如果我们的决定是保留(或拒绝)零假设，我们得出的结论是，不同组之间没有显著差异，或者变量之间没有关系。我们做出了正确的决定。如果零假设H0是否真的是假的，则组之间或变量之间存在显著差异，应予以拒绝。如果我们的决定是拒绝零假设H0≠？，我们得出的结论是，在组之间或变量之间存在显著的差异。我们又一次做出了正确的决定。然而，当您在不应该拒绝零假设的情况下拒绝该零假设时，就会发生类型I错误(概率=a)。类型II错误(概率=b)发生在您本应拒绝零假设却未能拒绝它的情况下。统计检验的功效定义为1-β，它是正确检测到总体差异的概率；或者检验的能力通常被定义为当零假设确实为假时拒绝零假设的概率。幂一般是可能分布的函数，通常由另一种假设下的参数确定。有几个因素影响功效。我们可以非正式地将这些因素分为定义权力因素和方法因素的参数。定义功率的参数以更“机械”的方式增加或减少功率。例如，对于两样本t检验，幂取决于总样本量、组样本量比率、α、平均差异或效应大小、标准差或变异性。幂、效果大小、样本大小和α是四个相互关联的参数，它们相互关联，使得每个参数都是其他三个参数的函数。换句话说，如果其中三个值是固定的，那么第四个值就完全确定了。因此，通过增加一个参数，可以减少(或增加)另一个参数：Alpha越高，功率越高(保持所有其他参数不变)；效果大小越大，需要的对象越少(给定相同的功率和Alpha级别)。尽管样本量计算可能会根据研究设计的类型而有所不同，但基本概念保持不变。
方法论因素包括实验设计、分组、统计程序和模型、时间点之间的相关性、反应变量和缺失数据等。方法论因素也通过更多的方法论问题来影响POWER。例如，重复测量设计实际上总是比横截面设计更强大；与较少的数据收集相比，响应变量的更多时间点也会增加功率。效应大小较小的两组通常需要比具有更大效果大小的那些更大的功率来检测。统计程序和模型也会影响功效：当模型的假设被违反时，非参数模型可能比参数模型更强大。相互作用项通常需要比主效应更强的能量来检测。有时，时间点之间相关性也会影响功率。响应变量的对数变换或对连续响应变量的归类也会导致能量的损失，因为调整模型会失去更多的能量，而Logistic或有序Logistic回归往往比普通的最小二乘回归需要更多的对象。最后，通常情况下，丢失数据会降低功耗，而不好的补偿方法会极大地降低功耗。在某些统计程序中，个案删除是处理丢失数据的默认设置。个案删除是失去权力的最大因素之一。因此，当您设计您的研究时，您需要尽一切可能将丢失的数据降至最低。
5.2 多样性差异使用T检验功效分析
①连续结果的幂公式:在微生物群研究中，样本中单个类群的数量和多样性可以用指数度量来概括，如Shannon多样性指数(均匀度)和Chao1指数(丰富度)。基本目标之一是比较不同群体的物种多样性。调查人员提出的问题是：一个群体是否比另一个群体存在更多的多样性。或者，来自两个数据集的alpha多样性的分析结果彼此一致？在统计假设框架下，可以提出问题。H0：组1中的分集与组2中的分集没有区别；Ha：组1的多样性与组2的不同。多样性度量可以形成假设检验、幂和样本量计算的基础。假设检验将在第7-12“章节”中介绍。重点介绍功效和样本量的计算。零假设是H0：u1=u2，其中u1和u2分别是组1和组2的真实分集平均值。另一种假设是Ha：u1不等于u2。零和可选的假设可以重写为H0：u1-u2=0和Ha：u1-u2= δ不等于0。功效分析的目标是以高概率检测大小δ的差异。一般功率公式由下式给出
对于连续端点，检验统计量为两组的平均差δ，然后
幂取决于α、δ、σ和n。如果σ已知，则在计算中使用正态分布，如果需要估计σ，则使用非中心t(或表)。假设σ已知，且n1=n2=n，那么我们可以使用数据来修改假设，如下所示：
如果出现以下情况，则拒绝H0分发
检验均数差异的一般样本量公式如下：
如果σ未知并且需要估计，并且假设n1=n2=n，那么我们可以将假设修改如下：如果出现以下情况，则拒绝H0分发
如果每组的样本量不平衡，n1不等于n2，我们可以推导出样本量公式如下
②ALS研究的多样性数据：为了说明使用t检验检验多样性差异的能力分析，我们比较了从肌萎缩侧索硬化症(ALS)研究中计算的G93A转基因小鼠的Shannon多样性，该研究是关于丁酸治疗对粪便微生物群的影响。G93A小鼠携带人类肌萎缩侧索硬化症(ALS)引起的SOD1突变，这些突变概括了肌萎缩侧索硬化症(ALS)患者的神经元和肌肉损伤。这些小鼠被广泛用于研究肌萎缩侧索硬化症的病理机制和试验治疗(张等人。2017年)。丁酸处理组9例，未处理组G93A突变体7例。计算了16只小鼠的香农多样性。图5.1中的分布显示，丁酸处理组产生更高的多样性，因为与不使用丁酸的对照组相比，该组的直方图向右移动(更高的多样性值)(注：垂直虚线标记每组的平均多样性)。在此数据中，丁酸处理组的香农多样性的平均值和标准差分别为2.504和0.170，而不加丁酸的对照组的香农多样性平均值和标准差分别为2.205和0.209。对于连续变量，效应大小正好是这两个均值的差值(δ=2.504-2.205)。
图5.1是通过以下R代码生成的。首先，我们加载分类群多度表，并将原始的按样本分类单元格式表转置为按分类单元分类单元格式表
> ##load abundance table
> abund_table=read.csv("ALSG93AGenus.csv",row.names=1,check.names=FALSE)
> abund_table_t<-t(abund_table)
然后，我们从vegan包中调用多样性函数来计算Shannon多样性。
> library(vegan)
> #use the diversity function (vegan package) to calculate Shannon index
> #make a data frame of Shannon index
> H<-diversity(abund_table_t, "shannon")
> df_H<-data.frame(sample=names(H),value=H,measure=rep("Shannon",length(H)))
第三，我们根据示例信息创建变量“Group”。
> #Obtain grouping information from sample data
> df_H$Group <- with(df_H,
+ ifelse(as.factor(sample)%in% c("A11-28F","A12-28F","A13-28F","A14-28F","A15-28F","A16-28F"),c("G93m1"),
+ ifelse(as.factor(sample)%in% c("A21-28F","A22-28F","A23-28F","A24-28F","A25-28F","A26-28F"),c("WTm1"),
+ ifelse(as.factor(sample)%in% c("C11-28F","C12-28F","C13-28F"),c("G93m4"),
+ ifelse(as.factor(sample)%in% c("C21-28F","C22-28F","C23-28F"),c("WTm4"),
+ ifelse(as.factor(sample)%in% c("B11-28F","B12-28F","B13-28F","B14-28F","B15-28F","D11-28F","D12-28F","D13-28F","D14-28F"),c("BUm3to3.5"),
+ c("NOBUm3to3.5")))))))
> df_H
sample value measure Group
A11-28F A11-28F 2.478 Shannon G93m1
A12-28F A12-28F 2.162 Shannon G93m1
A13-28F A13-28F 1.707 Shannon G93m1
A14-28F A14-28F 2.084 Shannon G93m1
A15-28F A15-28F 2.660 Shannon G93m1
A16-28F A16-28F 1.971 Shannon G93m1
A21-28F A21-28F 1.957 Shannon WTm1
...
整个数据集包括1个月、3个月、3.5个月和4个月的样本数据。我们感兴趣的是3到3.5个月期间治疗和不治疗的比较。因此，我们将3个月到3.5个月的数据子集如下
> library(dplyr)
> df_H_G6 <- select(df_H, Group,value)
> df_H_G93BUm3 <- filter(df_H_G6,Group=="BUm3to3.5"|Group=="NOBUm3to3.5")
> df_H_G93BUm3
Group value
1 BUm3to3.5 2.393
2 BUm3to3.5 2.677
3 BUm3to3.5 2.257
4 BUm3to3.5 2.700
5 BUm3to3.5 2.580
6 NOBUm3to3.5 2.262
7 NOBUm3to3.5 2.230
8 NOBUm3to3.5 2.433
9 NOBUm3to3.5 2.049
10 NOBUm3to3.5 1.814
11 BUm3to3.5 2.626
12 BUm3to3.5 2.334
13 BUm3to3.5 2.344
14 BUm3to3.5 2.625
15 NOBUm3to3.5 2.333
16 NOBUm3to3.5 2.312
最后，我们根据R代码生成图5.1。
> library(ggplot2)
> #split the plot into multiple panels
> p<-ggplot(df_H_G93BUm3, aes(x=value))+
+ geom_histogram(color="black", fill="black")+
+ facet_grid(Group ~ .)
> #Calculate the mean of each group
> ##calculate the average Shannon diversity of each group using the package Plyr
> library(plyr)
> mu <- ddply(df_H_G93BUm3, "Group", summarise, grp.mean=mean(value))
> head(mu)
Group grp.mean
1 BUm3to3.5 2.504
2 NOBUm3to3.5 2.205
> #add mean lines
> p+geom_vline(data=mu, aes(xintercept=grp.mean, color="red"),
+ linetype="dashed")
③使用R函数power.t.test()计算幂或样本量：为了检验香农多样性中没有差别的零假设，可以使用t检验或Wilcoxon秩和检验。我们把多样性的假设检验留到第7章。这里，我们重点说明如何使用R软件计算幂或样本量。在R中，可以使用Basic R中的函数power.t.test()和pwr包中的函数pwr.t.test()进行功率分析。我们使用power.t.test。此函数的用法如下所示：power.t.test (n = sample size, delta = effect size, sd = standard deviation, sig.level = 0.05, power = NULL, type = c(“two.sample”, “one.sample”, “paired”),alternative = (“two.sided”, “one.sided”))其中，n是每组样本大小的数目，delta 是均值的真差，sd是标准差，sig.level是显著性水平(类型I错误概率)，功率是测试的功率(1减去类型II错误概率)，类型是t测试的类型，并且备选是单侧或双侧测试。由于均值差的标准差未知，因此需要使用(5.5)中的公式进行估计。
> aggregate(formula = value * Group,
+ data = df_H_G93BUm3,
+ FUN = var)
Group value
1 BUm3to3.5 0.02892
2 NOBUm3to3.5 0.04349
> n 1 < - 9
> n2 <-7
> s1<-sqrt(0.02892)
> s2<-sqrt(0.04349)
> s=sqrt((n1-1)*s1^2+(n2-1)*s2^2)/(n1+n2-2)
> s
[1] 0.05012
现在可以通过调用power.t.test()函数获得功率。
> power.t.test(n=2:10,delta=2.504-2.205,sd=0.05012)
Two-sample t test power calculation
n = 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 1 0
delta = 0.299
sd = 0.05012
sig.level = 0.05
power = 0.8324, 0.9994, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000,
1.0000, 1.0000
alternative = two.sided
NOTE: n is number in *each* group
> df_P <-data.frame(n,power)
> df_P
n power
1 2 0.8324
2 3 0.9994
3 4 1.0000
4 5 1.0000
5 6 1.0000
6 7 1.0000
7 8 1.0000
8 9 1.0000
9 10 1.0000
从上面的功率分析可以看出，每组2只G93A小鼠的大小样本，随机分配到丁酸处理或不处理对照，将提供83%的功率来拒绝两组香农多样性没有差异的零假设。如果样本量增加到每组3个，则功率将增加到99%以上。我们可以生成功率和样本大小图，以可视化我们使用以下R代码拒绝零假设所需的功率和样本大小。
> n = c(2, 3, 4, 5, 6, 7, 8, 9, 10)
> power = c(0.8324, 0.9994, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000)
> power <- sapply(n, function (x) power.t.test(n=x, delta=2.504-2.205,sd=0.05012)$power)
> plot(n, power, xlab = "Sample Size per group", ylab = "Power to reject null",
+ main="Power curve for\n t-test with delta = 0.05",
+ lwd=2, col="red", type="l")
> abline(h = 0.90, col="blue")
实际上，幂和样本量的计算可以在几乎任何统计软件包中完成(例如，在R中，我们使用power.t.test(n=2：10，delta=2.504-2.205，sd=0.05012)来生成图5.2的数据)。假设比较组中的一个是预先指定的微生物群组，则需要单样本t检验来估计功率或样本大小。函数power.t.test()的默认测试是两个样本。您可以指定type=“one.sample“进行单样本功率分析，如下所示。
> power.t.test(n=2:10,delta=2.504-2.205,sd=0.05012, type = "one.sample")
5.3使用方差分析比较两个以上组的分集的功率分析
①单因素方差分析的假设和幂理论：对于涉及两个以上组的比较，可以使用单向方差分析(ANOVA)。例如，要比较三个或更多组中的多样性，假设可以是：H0：三个组的平均多样性相等；Ha：在这三个群体中，并不是所有的多样性都是平等的。差分析采用F检验。要在完全随机设计中计算一个处理或条件在k个水平上的全局F检验的功效，我们首先需要理解和定义假设。方差分析的基本思想是将多样性的总体方差划分为反映处理组或条件(因素水平)之间的差异和[由于测量误差(残差)]的处理组或条件内的差异的分量。对于出现在i1；…；K水平的因子a，每个水平有j1；…；J个观测值，典型的单因素方差分析模型可以表示为
根据统计假设框架，假设被定义为，
其中，ui=组i的平均值，k=组数，j=实验单位。均值相等的F检验采用单因素方差分析(ANOVA)假设数据是正态的，具有共同的组方差。也是N大于等于K，ni大于等于1，其中N是总样本量，ni是第i组的样本量。在零假设下，F统计量的分布服从中心F分布，而在另一种假设下，F统计量的分布服从带有非中心性参数k的非中心F分布。因此，当零假设为真时，它遵循中心F分布；当它为假时，它遵循非中心F分布。因此，幂可以定义为F统计量服从非中心分布的概率。 确切的功率是这样给出的，
上述功率用两个假设定义了两条分布曲线。零假设下的F分布定义了中心F分布，而在备选假设下的分布定义具有非中心性参数λ的相同F。R使用公式(5.10)计算幂。非中心性参数k是该公式中的关键参数。在R中λ是如何定义的？对于平衡设计，设n=平衡组的样本大小，则λ定义为
②使用R函数pwr.vova.test()计算幂或样本量：为了说明如何使用pwr.vova.test()获得不同样本大小的能力，我们使用以下过程。首先，我们将四组1个月和4个月的数据分成两组，分别为治疗组和对照组。
> df_H_G93WTm1N4 <- filter(df_H_G6,Group%in%c("G93m1","WTm1",
"G93m4","WTm4"))
> df_H_G93WTm1N4
Group value
1 G93m1 2.478
2 G93m1 2.162
3 G93m1 1.707
4 G93m1 2.084
5 G93m1 2.660
6 G93m1 1.971
7 WTm1 1.957
...
然后，通过拟合线性模型得到F统计量。
> fit = lm(formula = value*Group,data=df_H_G93WTm1N4)
> anova (fit)
Analysis of Variance Table
Response: value
Df Sum Sq Mean Sq F value Pr(>F)
Group 3 0.059 0.0195 0.23 0.88
Residuals 14 1.209 0.0863
最后，我们从pwr包中调用pwr.anova.test()函数来计算功率。
> library(pwr)
> pwr.anova.test(f= 0.23,k=4,n=45:55,sig.level=0.05)
Balanced one-way analysis of variance power calculation
k = 4
n = 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55
f = 0.23
sig.level = 0.05
power = 0.7276, 0.7383, 0.7486, 0.7586, 0.7683, 0.7777, 0.7868,
0.7956, 0.8041, 0.8123, 0.8202
NOTE: n is number in each group
上述结果表明，根据本次先导研究中检测到的效应大小，采用方差分析检验，每组需要52个样本才能获得80%的效果。
5.4 比较组间的分类单元
①比较比例的假设和基本幂和样本量公式：在微生物组数据集中，0表示样品中没有分类单元，1表示样品中存在分类单元。这样的二进制数据可以容易地转换成比例。如果研究对特定分类群的丰富度感兴趣，那么可以使用卡方检验来比较不同类群之间的比例。跨组比较单个指定分类单元的假设为 H0：组1中指定分类单元的百分比与组2中的没有区别；Ha：组1中指定分类单元的百分比与组2中的不同。方程式(5.4)可以很容易地修改以测试分类单元的比例。当我们测试分类单元比例时，我们正在评估对一个组(例如，treatment)响应的比例P1是否不同于对另一个组(例如， control）响应的比例P2。这等价于零假设H0：P1-P2=0与另一种假说Ha：P1-P2=δ≠0。
一般来说，当结果是二进制时，所需的样本大小如下所示。
类似地，在连续情况下，可以推导出两组比例的不平衡样本量公式。如果为n1≠n2，则两组样本大小之间的比率为r=n1÷n2。测试统计数据构造如下：
②使用R函数power.pro.test()进行功效分析：在我们的ALS研究中发现，口服2%丁酸钠2.5个月与丁酸产生菌Butyribrio fisolvens在物种水平上的丰度显著提高有关。我们有兴趣比较这个分类单元在处理组和对照组之间是否以不同的速率存在。为了设计未来的研究，我们希望计算样本量，以确保研究在此先导研究的基础上有足够的能力来检测效应大小。在本节中，我们将演示函数power.pro.test()，以便使用我们的ALS G93A小鼠数据进行功率分析，说明X2test和Fisher的精确测试。我们对丁状弧菌的丰度计数数据进行了转换将纤维隔离物转化为二元变量，以指示是否存在纤维隔离物丁状核弧菌。我们在表5.2中总结了转换后的数据。
下面的R代码加载丰度表
abund_table_Spe=read.csv("ALSG93AButyrivibrioSpecies.csv",row.names=1,check.names=FALSE)
原始丰度数据是一个以物种为行、以样本为列的计数表。下面的转置函数“t()”将该物种计数表转换成行的样本和列的物种。
> abund_table_Spe<-t(abund_table_Spe)
可以从样本标识符获得分组信息。
> grouping<-data.frame(row.names=rownames(abund_table_Spe),t(as.data.frame(strsplit(rownames(abund_table_Spe),"-"))))
将“B11”、“B12”、“B13”、“B14”、“B15”、“D11”、“D12”、“D13”、“D14”的9个样本随机分为丁酸治疗组，另外7个样本设为对照组。以下R代码用于创建一个组变量来表示分组信息。
> grouping$Group <-with(grouping,ifelse(as.factor(X1)%in%c("B11","B12","B13","B14","B15","D11","D12","D13","D14"),c("Butyrate"), c("Control")))
将物种丰度表和分组信息表合并成一个数据框后，即可进行数据分析。
> Butyrivibrio_G <-cbind(abund_table_Spe, grouping)
> rownames(Butyrivibrio_G)<-NULL
> Butyrivibrio_G
Butyrivibrio X1 X2 Group
1 14 B11 28F Butyrate
2 39 B12 28F Butyrate
3 18 B13 28F Butyrate
4 41 B14 28F Butyrate
5 19 B15 28F Butyrate
6 0 B21 28F Control
7 16 B22 28F Control
8 1 B23 28F Control
9 8 B24 28F Control
10 0 B25 28F Control
11 78 D11 28F Butyrate
12 4 D12 28F Butyrate
13 17 D13 28F Butyrate
14 94 D14 28F Butyrate
15 1 D21 28F Control
16 0 D22 28F Control
带0的丁状弧菌计数编码为“不存在”，否则编码为“存在”。
> Butyrivibrio_G$Present <- ifelse((Butyrivibrio_G$Butyrivibrio > 0), "Present","Absent")
> Butyrivibrio_G
Butyrivibrio X1 X2 Group Present
1 14 B11 28F Butyrate Present
2 39 B12 28F Butyrate Present
3 18 B13 28F Butyrate Present
4 41 B14 28F Butyrate Present
5 19 B15 28F Butyrate Present
6 0 B21 28F Control Absent
7 16 B22 28F Control Present
8 1 B23 28F Control Present
9 8 B24 28F Control Present
10 0 B25 28F Control Absent
11 78 D11 28F Butyrate Present
12 4 D12 28F Butyrate Present
13 17 D13 28F Butyrate Present
14 94 D14 28F Butyrate Present
15 1 D21 28F Control Present
16 0 D22 28F Control Absent
数据汇总为2×2列联表。
> library(MASS) # load the MASS package
> tbl = table(Butyrivibrio_G$Group, Butyrivibrio_G$Present)
> tbl # the contingency table
Absent Present
Butyrate 0 9
Control 3 4
上表显示了7个对照样品中有4个(57%)含有丁酸弧菌的比例，而9个丁酸盐样品中有9个(100%)含有丁酸治疗样本做到了。使用(5.16)和(5.17)中的上述公式，使用以下纯R码来实现样本大小和功效计算。
> p1=1.0
> p2=0.57
> r=1
> alpha=0.05
> beta=0.20
> (n2=(p1*(1-p1)/r+p2*(1-p2))*((qnorm(1-alpha/2)+qnorm(1-beta))/(p1-p2))^2)
[1] 10.4
> ceiling(n2)
[1] 11
> z=(p1-p2)/sqrt(p1*(1-p1)/n2/r+p2*(1-p2)/n2)
> (Power=pnorm(z-qnorm(1-alpha/2))+pnorm(-z-qnorm(1-alpha/2)))
[1] 0.8
不过，更方便的方法是使用下面的R函数power.pro.test。您可以指定多个样本来测试功效。
>power.prop.test(n=10:20, p1=1, p2=.57, sig.level=0.05, power=NULL, alternative=c("one.sided"), strict = FALSE)
Two-sample comparison of proportions power calculation
n = 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20
p1 = 1
p2 = 0.57
sig.level = 0.05
power = 0.7928, 0.8290, 0.8596, 0.8852, 0.9065, 0.9242, 0.9387, 0.9
506, 0.9603, 0.9682, 0.9746
alternative = one.sided
NOTE: n is number in *each* group
结果表明，在我们前期研究的基础上，每组11个样本可以获得83%的能量来检测效应大小。
③基于X2检验和Fisher精确检验的功率分析
比较比例的X2检验的功率理论：检验微生物群组成比例的假设也可以用卡方统计检验。X2test统计量被定义为，
与F统计量的分布类似，当零假设为真时，X2统计量的分布服从中心卡方分布。当它为假时，它遵循具有非中心性参数λ的非中心卡方分布。因此，本质上，卡方分布的幂是数据来自非中心卡方分布的概率。基本上，X2分布定义了两条曲线：一条是零假设下的中心X2分布，另一条是带有非中心性参数k的非中心X2分布，零假设下概率为0.05的卡方用垂直线表示。如果X2统计数据落在行的左侧，则表明它来自空分布，如果它在行的右侧，则它来自替代分布。线右侧的备选假设曲线下的面积表示测试的威力。例如，在微生物组研究中，我们假设组1或条件1的比例遵循中心X2分布，而组2或条件2遵循非中心X2分布。
Using the Function pwr.chisq.test()实现功效分析：为了找到指定功率所需的样本大小，R找到满足公式的N。在这一部分中，我们使用我们的丁状弧菌纤维隔离物存在/不存在的数据，使用来自ALS G93A小鼠数据集的数据，说明使用v2test和Fisher精确测试进行功率分析。要进行功率分析，我们首先需要估计影响的大小。对于连续测度，例如在香农分集的情况下，效应大小仅计算为均值的差值。在分类数据中，有几个度量可以用来作为效应大小，包括Cramer‘s V(φ系数u)、优势比和相对风险(Cohen，1988)。你可以根据你的学习兴趣选择一个。这里我们使用Cramer‘s V，它是名义变量关联性的度量。实际上，它是重新调整的皮尔逊卡方统计数据，其值在0到1之间。
其中，X2是零假设为真时的皮尔逊卡方，N是总样本大小，t是行数减一或列数减一中的较小者。如果r是行数，c是列数，则t=最小值(r−1，c−1)。当然，对于2X2列联表，这只是卡方除以观测次数的平方根，等于φ系数u。Cramer‘s V可以通过使用LSR软件包中的函数CramersV来估计。
> library(lsr)
> cramersV(tbl)
[1] 0.3833
给定Cramer‘s V，现在我们可以使用PWR包中的函数pwr.chisq.test()进行功率分析。
> library(pwr)
> pwr.chisq.test(w = 0.3833, N = 45:60, df = 1, sig.
level = 0.05, power = NULL)
Chi squared power calculation
w = 0.3833
N = 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60
df = 1
sig.level = 0.05
power = 0.7295, 0.7388, 0.7479, 0.7567, 0.7652, 0.7735, 0.7815,0.7893, 0.7969, 0.8042, 0.8113, 0.8182, 0.8248, 0.8313, 0.8375, 0.8435
NOTE: N is the number of observations
结果表明，每组需要54个样本才能在X2检验的基础上，以Cramer‘s V(φ系数u)为效应大小，以80%的权值正确拒绝零假设。
Using the Functions
power.fisher.test() and power.exact.test()进行功效分析：如果列联表中的单元格值较小(<5)，则卡方检验可能不正确，则应用Fisher精确检验。因此，我们将分别从statmod和Exact包中说明最近开发的两个R幂函数power.fisher.test()和power.exact.test()。与pwr.chisq.test()相比，两个版本的Fisher测试在一次测试中指定多个样本的灵活性较低。在下面的示例中，我们使用列联表中的相同比例(p1=1.0，p2=0.57)和每组15个样本来估计功率。
> library(statmod)
> power.fisher.test(1.0,0.57,15,15,alpha=0.05, nsim=1000)
[1] 0.844
以下R代码使用power.exact.test函数计算无条件精确测试的功率。
> library(Exact)
> power.exact.test(1.0, 0.57, 15, 15, method="Fisher")
$power
[1] 0.8454
$alternative
[1] "two.sided"
$method
[1] "fisher"
结果表明，每组需要15个样本才能在Fisher精确检验的基础上以84%的权值正确拒绝零假设。
5.5 利用Dirichlet-Polyomial模型比较所有类群的分布频率
①多元假设检验与Dirichlet-多项式模型：上面比较感兴趣的分类单元的方法是单变量的“一次一个分类单元”分析。由于微生物组数据具有多变量结构，我们经常通过跨组比较每个类群的丰度来进行多个类群的比较，然后进行调整以进行多次比较。单变量方法通常不如多变量方法强大，因为它没有考虑到分类群之间存在的相互作用，而多变量方法联合建模微生物类群丰度。微生物组研究人员一直在努力开发统计方法，以进行关于处理或实验因素对整个细菌分类群组合的影响的多变量假设，并估计此类实验的样本量。在比较微生物群的多变量假设检验框架下，无效假设之一可以是：H0：不同的处理或实验因素对整个类群的影响没有不同。Ha：不同的处理或实验因素对整个类群的影响是不同的。多变量统计方法之一是Dirichlet多项式分布，它已被证明能很好地模拟微生物组数据。La Rosa及其同事对Dirichlet多项式模型进行了重新参数化，使其适用于基于位置(均值比较)和尺度(方差比较/离散度)之间的差异进行跨组假设检验。该方法具有参数估计、多变量假设检验、幂和样本量计算等功能。在这一章中，我们首先介绍了重新参数化的Dirichlet多项式模型，然后重点介绍了使用HMP软件包计算幂和样本量。
这种结构的计数数据一般采用多项式分布进行分析。然而，它使用多项分布的适宜性假设每个类别的真实频率在所有样本中都是相同的(La Rosa等人。(2012年)。此外，当数据呈现过度离散性时，多项式模型可能会导致I类误差增加。因此，使用多项式分布来模拟微生物组计数数据是不合适的，因为由于样本的可变性，微生物组数据中的每个分类单元在所有样本中并不相同。为了将微生物群分类单元计数数据与该数据结构相匹配，重新参数化了Dirichlet-多项分布。其特征在于两组参数：π和⊙：
其中，⊙是超色散参数，或色散的量度。它测量样本内相对于多项分布的变异性过剩；⊙≥0表示过度分散。与多项分布相比，Dirichlet多项分布的参数化使得它适合于通过比较π向量和比较⊙值来进行基于位置差异的跨组假设检验。例如，对于两组比较(即对照和病例)，零假设：π1=π2表示群落组成相等。有人建议使用Koehler和Wilson广义Wald检验来检验零假设。作为Dirichlet-Polyomial，重新参数化的Dirichlet-Polyomial模型既可用于分析分类单元组成数据，也可用于分析秩丰度分布(RAD)数据。这两种方法分别称为“分类群组成数据分析”和“秩丰度分布数据分析。两种方法有不同的侧重点：一种是群落组成(那里有什么细菌)；另一种是群落结构(如丰富度和多样性)。
②Dirichlet-多项式模型下的幂和样本量计算：Dirichlet多项式模型的幂和样本量取决于概率模型参数、假设正在测试中，以及效果的大小。模拟研究表明，效应大小、超分散和样本量影响功率。功率随效应尺寸的增大或样本量的增大而增大，随过色散的减小而增大。在一些示例中，在保持效应大小、超色散和样本大小恒定的情况下，读取次数也影响功率，功率随着读取次数的增加而增加。总之，统计能力受到每个处理的样本数量(样本大小或复制水平)和每个样本的序列读取数量(测序深度)的影响，以及过度分散。为简单起见，让我们考虑一下在两个实验组之间比较分类群频率的微生物群样本的问题，例如用丁酸处理的小鼠和不用丁酸处理的小鼠。为了将分类群频率与先前指定的微生物种群、多于两个组以及其他假设检验进行比较。对于两个组之间的π的比较，假设每个组的模型参数π和⊙都是已知的，则检验假设形成为：H0: π1=π2 versus the alternative HA: π1≠π2
影响的大小由分类频率π1和π2的矢量彼此之间的距离来定义。有几种方法可以测量效果大小。其中之一是(5.20)中的Cramer‘s V。对于表5.3中的微生物组数据结构，行r是要比较的组的数量，列c是分类群的总数。根据公式，Cramer‘s V取决于样本大小和读取次数。为了减少样本大小和读取次数的影响，提出了一个修正的Cramer u准则，范围从0到1，前者表示两个类群的分类单元频率相同，后者表示分类单元频率差异最大。
③使用HMP软件包计算功率和尺寸：HMP是一个R包，用于假设检验和能力计算，用于比较最初来自HMP的元基因组样本(La Rosa等人。2016)。功率和样本量的计算可以通过使用统计量的渐近分布或通过HMP软件包实现的Monte-Carlo模拟来进行。在这里，我们使用ALS G93A小鼠的例子，通过蒙特卡罗模拟说明功率和样本大小的计算。
准备数据集以使用HMP包：首先，安装并调用HMP包
> install.packages("HMP",repo="http://cran.r-project.org", dep=TRUE)
> library(HMP)
然后设置工作目录，加载类群丰度数据集
> setwd("F:/Home/MicrobiomeStatR/Analysis")
> Buty=read.csv("ALSG93A3.5mButyrateGenus.csv",row.names=1,check.names=FALSE)
> NOButy=read.csv("ALSG93A3.5mNoButyrateGenus.csv",row.names=1,check.names=FALSE)
原始数据集具有按样本分类的格式。如表5.3所示，重新参数化的Dirichlet多项式模型和HMP包需要数据集具有分类抽样格式。因此，我们使用t()函数将原始数据集转置为具有Sample by Taxa格式的数据集。
> head(Buty)
B11-28F B12-28F B13-28F B14-28F B15-28F
Lactococcus 17 204 8 7 4
Tannerella 646 170 670 421 548
Barnesiella 6 2 12 8 21
Bacteroides 604 406 436 260 443
Hydrogenoanaerobacterium 1 0 7 1 0
Clostridium 179 398 564 400 737
> head(NOButy)
B21-28F B22-28F B23-28F B24-28F B25-28F
Lactococcus 52 21 0 9 1
Tannerella 787 756 395 1266 1111
Barnesiella 12 21 8 24 17
Bacteroides 130 241 192 228 315
Hydrogenoanaerobacterium 0 0 0 0 1
Clostridium 458 344 418 167 334
> Buty_t <- t(Buty)
> NOButy_t<-t(NOButy)
> head(Buty_t)
Lactococcus Tannerella Barnesiella Bacteroides
B11-28F 17 646 6 604
B12-28F 204 170 2 406
B13-28F 8 670 12 436
B14-28F 7 421 8 260
B15-28F 4 548 21 443
> head(NOButy_t)
Lactococcus Tannerella Barnesiella Bacteroides
B21-28F 52 787 12 130
B22-28F 21 756 21 241
B23-28F 0 395 8 192
B24-28F 9 1266 24 228
B25-28F 1 1111 17 315
以下代码检查每个数据集中的分类群数和样本数
> ncol(Buty_t) # for the number of taxa
[1] 196
> nrow(Buty_t) # for the number of samples
[1] 5
> ncol(NOButy_t) # for the number of taxa
[1] 196
> nrow(NOButy_t) # for the number of samples
[1] 5
利用分类群组成数据分析计算功率和大小：函数MC.Xdc.Statistics()使用模拟和似然比检验来提供几个样本Dirichlet-多项式参数检验比较的幂和大小。这里，我们比较“Buty”和“NOButy”的两个样本集。一个示例语法如下：
MC.Xdc.statistics(group.Nrs, numMC = 1000, alphap, type = “ha”, siglev =0.05, est = “mom”)
其中，group.Nrs用于指定组中每个样本的读取次数/序列深度。NumMC用于指定蒙特卡罗实验的次数。在实践中，至少应该指定1000个。如果alphap用于计算测试统计数据的大小(类型I错误)，即当指定type=“hnull”时，则alphap指定矩阵，其中行是参考组的alpha参数的向量。例如，以下代码：Alphap<Fit_Buty$Gamma指定参考组的“Buty”伽马矩阵。如果alphap用于测试统计数据的计算能力(类型II错误)，即默认值，或者当指定type=“ha”时，则alphap表示由每个组中每个分类单元的alpha参数矢量组成的矩阵。例如，以下代码：alphap<-rbind(Fit_Buty$Gamma，Fit_NOButy$Gamma)为每个组中的每个分类群指定“Buty”和“NOButy”伽马矩阵。如果type=“hnull”，则计算测试统计信息的大小；如果type=“ha”，则计算测试统计信息的幂，这也是默认设置。Siglev用于指定测试或功率计算大小的重要级别。默认值为0.05。est用于指定要与似然比检验统计量“mle”或“mom”一起使用的参数估计器的类型。默认值为“mom”。HMP软件包的作者指出，“最大似然估计”的运行时间要长得多，并且在小样本情况下不是最优的，而“mom”的结果在小样本情况下更为保守。首先，使用函数DM.MoM()获取数据的Dirichlet多项式参数列表(即loglik、Gamma、pi和theta)。
> fit_Buty <- DM.MoM(Buty_t);fit_NOButy <- DM.MoM(NOButy_t)
> fit_Buty
其次，设置蒙特卡罗实验的次数，这里我们使用1000，这是实践中推荐的最小值，并生成每个样本的读取次数。
> numMC <- 1000
> ##The first number is the number of reads and the second is the number of samples or subjects
> nrsGrp1 <- rep(1000, 10)
> nrsGrp2 <- rep(1000, 10)
> group_Nrs <- list(nrsGrp1, nrsGrp2)
第三，计算测试统计数据的大小(I类错误)。
> alphap <- fit_Buty $gamma
> pval1 <- MC.Xdc.statistics(group_Nrs, numMC, alphap, ”hnull”)
> pval1
[1] 0.000999
最后，测试统计数据的计算能力(类型II错误)。
> alphap <- rbind(fit_Buty$gamma, fit_NOButy$gamma)
> pval2 <- MC.Xdc.statistics(group_Nrs, numMC, alphap)
> pval2
[1] 0.999
表5.3显示了基于从Buty和NOButy 10样本数据集中获得的Dirichlet多项式参数，使用5%显著性水平比较Buty和NOButy种群分类群频率的功率分析。每个条目表示指定采样数和0.05显著性水平下的读取数所获得的功率。例如，对于样本数=10，并且每个样本的读取次数=500，研究具有98.8%的能力来检测在数据中观察到的效应大小。结果显示，当样本较小时，增加读取次数(在本例中为5)会影响功率；但当达到较大样本大小时(在此情况下，为10)，增加读取次数不会影响功率。结果还表明，当获得足够的样本时，增加样本量并不会增加功率。
利用秩丰度分布数据分析计算功率和尺寸：几种样本RAD-概率均值检验与已知比例参考向量的比较
函数MC.Xmc.Statistics()使用模拟和广义Wald-type统计量来提供几个样本RAD概率均值测试与已知比例参考向量比较的能力和大小。这里，我们再次比较“Buty”和“NOButy”的两个样例集。一个示例语法如下所示：MC.Xmc.statistics(group.Nrs, numMC = 1000, pi0, group.pi, group.theta, type =“ha”, siglev = 0.05)。其中，group.Nrs、numMC、type和siglev如MC.Xdc.Statistics()函数中定义。Pi0是RAD概率平均向量。如果group.pi=“hnull”，则忽略此参数；如果group.pi=“ha”，则指定一个矩阵，其中每行都是每个组的向量pi值。Group.theta是每个组的超离散值的向量。分析稀有类群的技术挑战是汇聚率低，测试统计不精确。为了提高收敛速度和准确估计，在微生物群研究中，一种方法是去掉频率较低的分类群，例如在任何百分比小于1的样本中；另一种方法是将所有频率较低的分类群合并起来。分类单元，例如，将低于1%的所有类群加权平均为一个分类单元，称为“其他”分类单元。在上述分类群组成数据分析中，我们在调用MC.Xdc.Statistics()函数时没有使用这些技术。在这里，我们使用Data.filter()函数按丰度递减的顺序对分类单元进行排序，并将不太丰富的分类单元折叠成一个标记为“Other”的类别。一个示例语法如下所示：Data.filter(data, order.type = “sample”, minReads = 1000, numTaxa = 10,perTaxa = NULL)。其中，data是每个样本(行)的分类计数(列)矩阵。如果order.type=“sample”，则根据分类频率对分类进行排名；如果order.type=“data”，则根据所有样本的累计分类计数对分类进行排名(默认)。假设minReads=一个读取截止值，则读取总数小于该读取截止值的样本将被删除。参数numTaxa和perTaxa，一次调用中只能指定一个。参数numTaxa用于指定要保留的分类群数量，同时将其他(不太丰富的)分类群折叠为“其他”类别。参数perTaxa用于合并要保留的数据的百分比，同时折叠剩余的分类群。首先，使用Data.filter()函数按丰度递减的顺序对分类群进行排序，并将丰度较低的分类群折叠为“其他”类别。
> filter_Buty<- Data.filter(Buty_t, "sample", 1000, 10)
> head(filter_Buty)
Other
B11-28F 646 604 265 196 179 159 103 53 52 32 261
B12-28F 406 398 312 242 204 170 102 82 50 45 382
B13-28F 670 564 436 110 65 62 59 58 47 39 266
B14-28F 421 400 260 149 111 67 64 58 49 48 421
B15-28F 737 548 443 281 214 97 94 69 59 53 476
> filter_NOButy<- Data.filter(NOButy_t, "sample", 1000, 10)
> head(filter_NOButy)
Other
B21-28F 787 458 231 130 114 110 67 61 52 35 220
B22-28F 973 756 344 312 241 151 145 56 40 38 256
B23-28F 418 395 192 165 80 42 41 37 36 34 238
B24-28F 1266 425 402 228 167 142 113 54 42 25 178
B25-28F 1111 334 315 102 59 53 43 35 22 19 133
其次，使用函数DM.MoM()获取数据的Dirichlet多项式参数列表(即loglik、Gamma、pi和theta)。
> fit_Buty <- DM.MoM(Buty_t);fit_NOButy <- DM.MoM(NOButy_t)
> fit_Buty$pi
0.23155 0.20212 0.13796 0.07863 0.06215 0.04462 0.03393 0.02573
Other
0.02066 0.01745 0.14520
> fit_NOButy$pi
0.36373 0.18909 0.11850 0.07482 0.05278 0.03977 0.03266 0.01940
Other
0.01533 0.01206 0.08185
> fit_Buty$theta
[1] 0.007523
> fit_NOButy$theta
[1] 0.01615
第三，设置蒙特卡罗实验的次数，这里我们使用1000，这是实践中推荐的最小值，并生成每个样本的读取次数。
> numMC <- 1000
> #The first number is the number of reads and
> #the second is the number of subjects
> nrsGrp1 <- rep(1000, 10);nrsGrp2 <- rep(1000, 10)
> group_Nrs <- list(nrsGrp1, nrsGrp2)
第四，设置各类群的类群频率向量(类群比例)和超分散参数的值。
> pi0 <- fit_Buty$pi
> group_theta <- c(0.007523, 0.01615)
第五，计算测试统计数据的大小(I类错误)。
> pval1 <- MC.Xmc.statistics(group_Nrs, numMC, pi0, group.theta=group_theta, type="hnull")
> pval1
[1] 0.08492
最后，计算测试统计数据的功率(类型II错误)。
> group_pi <- rbind(fit_Buty$pi, fit_NOButy$pi)
> pval2 <- MC.Xmc.statistics(group_Nrs, numMC, pi0, group_pi, group_theta)
> pval2
[1] 0.999
表5.4显示了使用函数MC.Xmc.Statistics()进行秩丰富分布(RAD)数据分析的功率分析结果。与分类群组成数据分析相比，RAD数据分析获得了更高的检测效果大小的能力。例如，每组仅5个样本，每个样本500个读取，我们可以达到89.61%的功率。实际上，与分类群组成数据分析相比，RAD数据分析方法降低了功耗，因为当忽略分类群标签时，数据中的信息会丢失。在这种情况下获得的更高的能力是因为RAD数据分析将196个分类群减少到11个分类群，其中包括1个“不那么频繁地汇集在一起的丰富的分类群”。
几种样本比例向量未知的RAD-概率均值检验比较：函数MC.Xmcupo.Statistics()使用模拟和广义的Wald-type统计量来提供几个样本RAD概率均值测试比较的功率和大小，而不需要比例参考向量。这里，我们再次比较“Buty”和“NOButy”的两个样例集。一个示例语法如下所示：MC.Xmcupo.statistics(group.Nrs, numMC = 1000, pi0, group.pi, group.theta,type = “ha”, siglev = 0.05)。其中，group.Nrs、numMC、type、siglev、pi0、group.pi和group.theta在函数MC.Xmc.Statistics()中定义。函数MC.Xmcupo.Statistics()的运行方式与函数MC相同。Xmc.Statistics()。
> ##Generate the number of reads per sample
> ##The first number is the number of reads and the second is the number of Samples or subjects
> nrsGrp1 <- rep(1000, 10);nrsGrp2 <- rep(1000, 10)
> group_Nrs <- list(nrsGrp1, nrsGrp2)
> pi0 <- fit_Buty$pi
> group_theta <- c(0.007523, 0.01615)
> ##Computing size of the test statistics (Type I error)
> group_theta <- c(fit_Buty$theta, fit_NOButy$theta)
> pval1 <- MC.Xmcupo.statistics (group_Nrs, numMC, pi0, group.theta=group_theta, type="hnull")
> pval1
[1] 0.004995
> ##Computing Power of the test statistics (Type II error)
> group_pi <- rbind(fit_Buty$pi, fit_NOButy$pi)
> pval2 <- MC.Xmcupo.statistics(group_Nrs, numMC, group.pi=-group_pi, group.theta=group_theta)
> pval2
[1] 0.9231
表5.5显示了使用函数MC.Xmcupo.Statistics()进行秩丰富分布(RAD)数据分析的功率分析结果。函数MC.Xmcupo.Statistics()实现的幂小于函数MC.Xmc.Statistics()的幂，这是因为没有使用比例的参考向量。
④使用HMP程序包计算效果：HMP包具有计算效果大小的功能。这里，我们使用函数Xmcupo.ffectsize()来说明效果大小的计算，该函数计算测试统计量Xmcupo.sevsample()的Cramer‘s Phi和修改后的Cramer’s Phi Criteria。使用函数Xmcupo.sevsample()进行假设测试。在这里，我们将重点放在效果大小的计算上。语法为：Xmcupo.effectsize(group.data)其中，group.data是一个列表，其中每个元素都是每个样本(行)的分类计数(列)矩阵。
> ##Combine the data sets into a single list
> group_data <- list(filter_Buty, filter_NOButy)
> effect <- Xmcupo.effectsize(group_data)
> effect
Chi-Squared Cramer Phi Modified-Cramer Phi
20.97915 0.02899 0.15208
P value
0.02124
在基于检验统计量Xmcupo.sev样本的RAD数据分析中观察到的效应大小分别为：Cramerφ0.03，Modified-Cramerφ0.15。
5.6 总结
本章重点介绍微生物组数据的功率和样本量计算。介绍了单变量和多变量分析方法。对于单变量分析方法，我们介绍并说明了参数和非参数检验和模型；对于多变量分析方法，重点介绍了使用参数Dirichlet多项式模型计算幂和样本量。在每个分析中都引入了统计理论和相关的R软件包。我们从介绍假设检验和力量分析的要素开始，并利用微生物组数据形成假设检验和力量分析。其次，我们划分了四个部分，分别涵盖力量分析和样本量分析，分别使用t-test检验多样性差异；使用方差分析比较两个以上组的多样性；比较组间感兴趣的分类群；以及使用Dirichlet-Polyomial模型比较组间所有类群的频率。微生物组能力分析对于微生物组研究的研究设计是非常重要的。文献中的假设检验方法和模型有限，特别是针对合适的多变量模型的假设检验更是凤毛麟角。
