5.1假設(shè)檢驗(yàn)和功效分析
①假設(shè)檢驗(yàn)：統(tǒng)計(jì)學(xué)的主要目的是根據(jù)樣本信息對未知的總體參數(shù)進(jìn)行推斷：從觀察到的樣本數(shù)據(jù)中得出關(guān)于總體的結(jié)論赦拘。統(tǒng)計(jì)推斷包括總體參數(shù)的估計(jì)和關(guān)于總體參數(shù)的假設(shè)檢驗(yàn)。假設(shè)檢驗(yàn)是一種用來檢驗(yàn)統(tǒng)計(jì)假設(shè)的統(tǒng)計(jì)程序玄叠。這是制定統(tǒng)計(jì)程序以決定接受或拒絕假設(shè)的統(tǒng)計(jì)推理的主要領(lǐng)域序愚，而冪和樣本量的計(jì)算與假設(shè)檢驗(yàn)有關(guān)率碾。統(tǒng)計(jì)假設(shè)是關(guān)于人口的斷言或理論，特別是關(guān)于人口分布的參數(shù)帆吻，如位置(均值)域那、尺度(方差/離散度)。本章的目的是在一般情況下，特別是在微生物群研究中次员，介紹假設(shè)檢驗(yàn)和相關(guān)的功率和樣本量計(jì)算败许。
統(tǒng)計(jì)測試的要素：統(tǒng)計(jì)檢驗(yàn)有四個(gè)要素：(I)零假設(shè)，H0淑蔚；(Ii)替代假設(shè)市殷，H1；(Iii)檢驗(yàn)統(tǒng)計(jì)刹衫，和(Iv)拒絕區(qū)域醋寝。每個(gè)假設(shè)檢驗(yàn)程序由兩個(gè)假設(shè)組成。第一個(gè)假設(shè)是零假設(shè)(用H0表示)带迟，這是一種關(guān)于一個(gè)或多個(gè)總體參數(shù)的特定值的理論音羞。理論通常表述為H0：參數(shù)=特定值，例如H0：l仓犬？0：25嗅绰。我們想要根據(jù)樣本中包含的信息來檢驗(yàn)的假設(shè)，通常是在希望它被拒絕的情況下制定的搀继。第二個(gè)假設(shè)被稱為替代假設(shè)(用Ha表示)窘面，這是一種通常以與零假設(shè)相反的形式提出的理論。通常叽躯，零假設(shè)聲稱沒有差異或沒有變化财边，而另一種假設(shè)則聲明有一些差異或變化已經(jīng)發(fā)生。因此点骑，零假設(shè)總是聲明總體參數(shù)與索賠值完全相等酣难，而另一種假設(shè)允許該參數(shù)有幾個(gè)不同的值。例如黑滴，在一項(xiàng)臨床試驗(yàn)中鲸鹦，如果研究人員想知道與接受安慰劑20 mm的患者相比，服用新藥的患者是否降低了平均標(biāo)準(zhǔn)化平均得分的0.20分以下跷跪。零假設(shè)將是H0：l？0：20齐板，而替代假設(shè)將是以下之一：HA：l[0：20吵瞻，其中l(wèi)表示總體的平均標(biāo)準(zhǔn)化平均分?jǐn)?shù)。檢驗(yàn)統(tǒng)計(jì)量是用于決定是否拒絕零假設(shè)的樣本統(tǒng)計(jì)量甘磨。它是樣本觀測值和一些已知常數(shù)的函數(shù)橡羞，統(tǒng)計(jì)決策將基于這些常數(shù)。拒絕區(qū)域指定將拒絕零假設(shè)的測試統(tǒng)計(jì)量的數(shù)值济舆。具體地說卿泽，它是一個(gè)概率區(qū)域，它告訴我們我們的理論或假設(shè)可能是真的，也可能不是真的签夭。那么問題是：如何選擇拒絕區(qū)域齐邦？簡而言之，排斥區(qū)域與α水平相關(guān)第租。您需要選擇您愿意接受的Alpha級別措拇。例如，如果您想要95%確信您的結(jié)果是重要的慎宾，您將選擇5%的alpha級別(100%-95%)丐吓。“5%水平”是拒絕區(qū)域趟据。
檢驗(yàn)統(tǒng)計(jì)假設(shè)的步驟：從樣本統(tǒng)計(jì)到總體參數(shù)估計(jì)券犁，統(tǒng)計(jì)假設(shè)檢驗(yàn)起著基礎(chǔ)性的作用。您可以使用兩種方法進(jìn)行假設(shè)檢驗(yàn)：使用臨界值和使用p值汹碱。對于臨界值粘衬，推理步驟如下：(1)由于初始研究假設(shè)未知，說明相關(guān)的無效假設(shè)和備選假設(shè)比被；(2)對測試樣本進(jìn)行統(tǒng)計(jì)假設(shè)色难，如統(tǒng)計(jì)獨(dú)立性或觀測值的分布；(3)選擇合適的檢驗(yàn)等缀，并指定相應(yīng)的檢驗(yàn)統(tǒng)計(jì)量T枷莉；(4)從假設(shè)出發(fā)，推導(dǎo)出零假設(shè)下檢驗(yàn)統(tǒng)計(jì)量的分布尺迂；(5)選擇通常為小值(例如0.01笤妙、0.05)的顯著性水平(臨界區(qū))a，其被稱為檢驗(yàn)的顯著性水平或低于其將拒絕零假設(shè)的概率閾值噪裕；(6)由觀測值計(jì)算檢驗(yàn)統(tǒng)計(jì)量T的觀測值蹲盘；(7)決定要么拒絕零假設(shè)，轉(zhuǎn)而支持備選方案膳音，要么接受它召衔。您還可以決定拒絕或接受具有p值的零假設(shè)。如果要使用p值方法祭陷，則需要計(jì)算p值苍凛。這是在零假設(shè)下，獲得等于或比實(shí)際觀察到的結(jié)果更極端的結(jié)果的概率兵志。如果p值落在拒絕區(qū)域醇蝴，則意味著您具有統(tǒng)計(jì)上顯著的結(jié)果；因此想罕，您可以拒絕零假設(shè)悠栓，并得出結(jié)論，認(rèn)為替代假設(shè)為真。如果p值落在拒絕區(qū)域之外惭适，則意味著您的結(jié)果不能提供足夠的證據(jù)來反對無效假設(shè)笙瑟。
微生物組數(shù)據(jù)中的假設(shè)檢驗(yàn)：對于微生物群落比較，一般假設(shè)測試可以寫成：(1) H0:?π1=π0 versus HA:?π1不等于π0
(2) H0:?π1=...πj=...πJ versus HA:?πi不等于πj腥沽。上述兩種假設(shè)檢驗(yàn)分別類似于經(jīng)典統(tǒng)計(jì)學(xué)中的單樣本t檢驗(yàn)和兩樣本t檢驗(yàn)或方差分析逮走，它們構(gòu)成了微生物群研究中用于比較的基本統(tǒng)計(jì)假設(shè)檢驗(yàn)框架。A.分類群與先前指定的微生物群落的平均比例 B.兩個(gè)樣本集分類群的平均比例 C.來自兩個(gè)以上類群的平均分類群比例 D.分類群的平均比例和比例都是跨組的今阳。
微生物組研究中的功率和樣本量計(jì)算基于假設(shè)檢驗(yàn)框架师溅。微生物組數(shù)據(jù)既有其獨(dú)特性，又有其他研究領(lǐng)域的共性盾舌。因此墓臭，假設(shè)檢驗(yàn)和冪以及樣本量計(jì)算具有相似的設(shè)置。對于微生物組數(shù)據(jù)的共同特征妖谴，根據(jù)這些數(shù)據(jù)值的分布方式和要比較的組數(shù)窿锉，可以對微生物組假設(shè)使用標(biāo)準(zhǔn)t檢驗(yàn)、方差分析(ANOVA)或相應(yīng)的非參數(shù)檢驗(yàn)膝舅。對于微生物組數(shù)據(jù)的獨(dú)特特征嗡载，研究人員試圖開發(fā)適當(dāng)?shù)慕y(tǒng)計(jì)分析工具，包括功率和大小計(jì)算仍稀，以更好地?cái)M合數(shù)據(jù)洼滚。例如，考慮到微生物組數(shù)據(jù)結(jié)構(gòu)技潘，開發(fā)了多變量分析工具來考慮分類群之間的相互作用或相互關(guān)系遥巴。其中，HMP軟件包采用基于分類群計(jì)數(shù)的Dirichlet-多項(xiàng)式模型的參數(shù)方法享幽。
②功效分析與樣本量計(jì)算
功效分析與樣本量計(jì)算的重要性：功率分析最常見的目的是確定合理檢測給定大小的影響所需的最小對象铲掐。此外，如果有可用的效果大小和受試者數(shù)量值桩，則使用功率分析來確定功率摆霉。功率分析還可以用來比較不同的統(tǒng)計(jì)檢驗(yàn)過程，例如奔坟，在同一假設(shè)的參數(shù)檢驗(yàn)和非參數(shù)檢驗(yàn)之間斯入。例如，當(dāng)你為NIH撥款提案計(jì)劃一項(xiàng)研究或設(shè)計(jì)一項(xiàng)實(shí)驗(yàn)時(shí)蛀蜜，你會問一個(gè)常見的問題：“我們需要多少科目？”這個(gè)問題很重要增蹭，因?yàn)闃颖敬笮?yīng)該足夠大滴某，如果存在科學(xué)興趣的效應(yīng)大小，它就有很好的機(jī)會被檢測出來。以小樣本進(jìn)行研究霎奢，在時(shí)間和金錢上都是浪費(fèi)資源户誓，因?yàn)榻Y(jié)果往往是沒有定論的。但是幕侠，也不建議使用大樣本量帝美。首先，它不會有什么用處晤硕，因?yàn)樵诳茖W(xué)上重要性不大的效應(yīng)大小在統(tǒng)計(jì)上也可能是有意義的悼潭。其次，出于經(jīng)濟(jì)原因舞箍，確定合適的樣本量是很重要的舰褪。如果能從較小的樣本中準(zhǔn)確地找到答案，那就是浪費(fèi)有限的可用資源疏橄。第三占拍，在人體研究或隨機(jī)對照試驗(yàn)中，招募比所需更多的受試者可能是不道德的捎迫。參與治療的患者不應(yīng)被誤用晃酒，特別是安慰劑組患者。
功效分析：功效是當(dāng)零假設(shè)(H0)實(shí)際上為假或替代假設(shè)(HA)為真時(shí)窄绒，測試正確拒絕該零假設(shè)(H0)的概率贝次。它可以等效地被認(rèn)為是當(dāng)備選假設(shè)(HA)為真時(shí)接受該假設(shè)的概率；也就是說颗祝，能力是測試檢測效果的能力浊闪，如果該效果確實(shí)存在的話。因此螺戳，要理解權(quán)力的概念搁宾，我們必須理解假設(shè)檢驗(yàn)的原理。用于統(tǒng)計(jì)假設(shè)檢驗(yàn)的四種可能的決策結(jié)果及其相關(guān)概率之間的關(guān)系如表5.1所示倔幼。
如果零假設(shè)H0是否真的為真盖腿，則組間無顯著差異或變量間無關(guān)系，不應(yīng)予以否定损同。如果我們的決定是保留(或拒絕)零假設(shè)翩腐，我們得出的結(jié)論是，不同組之間沒有顯著差異膏燃，或者變量之間沒有關(guān)系茂卦。我們做出了正確的決定。如果零假設(shè)H0是否真的是假的组哩，則組之間或變量之間存在顯著差異等龙，應(yīng)予以拒絕处渣。如果我們的決定是拒絕零假設(shè)H0≠？蛛砰，我們得出的結(jié)論是罐栈，在組之間或變量之間存在顯著的差異。我們又一次做出了正確的決定泥畅。然而荠诬，當(dāng)您在不應(yīng)該拒絕零假設(shè)的情況下拒絕該零假設(shè)時(shí)，就會發(fā)生類型I錯(cuò)誤(概率=a)位仁。類型II錯(cuò)誤(概率=b)發(fā)生在您本應(yīng)拒絕零假設(shè)卻未能拒絕它的情況下柑贞。統(tǒng)計(jì)檢驗(yàn)的功效定義為1-β，它是正確檢測到總體差異的概率障癌；或者檢驗(yàn)的能力通常被定義為當(dāng)零假設(shè)確實(shí)為假時(shí)拒絕零假設(shè)的概率凌外。冪一般是可能分布的函數(shù)，通常由另一種假設(shè)下的參數(shù)確定涛浙。有幾個(gè)因素影響功效康辑。我們可以非正式地將這些因素分為定義權(quán)力因素和方法因素的參數(shù)。定義功率的參數(shù)以更“機(jī)械”的方式增加或減少功率轿亮。例如疮薇，對于兩樣本t檢驗(yàn)，冪取決于總樣本量我注、組樣本量比率按咒、α、平均差異或效應(yīng)大小但骨、標(biāo)準(zhǔn)差或變異性励七。冪、效果大小奔缠、樣本大小和α是四個(gè)相互關(guān)聯(lián)的參數(shù)掠抬，它們相互關(guān)聯(lián)，使得每個(gè)參數(shù)都是其他三個(gè)參數(shù)的函數(shù)校哎。換句話說两波，如果其中三個(gè)值是固定的，那么第四個(gè)值就完全確定了闷哆。因此腰奋，通過增加一個(gè)參數(shù)，可以減少(或增加)另一個(gè)參數(shù)：Alpha越高抱怔，功率越高(保持所有其他參數(shù)不變)劣坊；效果大小越大，需要的對象越少(給定相同的功率和Alpha級別)屈留。盡管樣本量計(jì)算可能會根據(jù)研究設(shè)計(jì)的類型而有所不同局冰，但基本概念保持不變括儒。
方法論因素包括實(shí)驗(yàn)設(shè)計(jì)、分組锐想、統(tǒng)計(jì)程序和模型、時(shí)間點(diǎn)之間的相關(guān)性乍狐、反應(yīng)變量和缺失數(shù)據(jù)等赠摇。方法論因素也通過更多的方法論問題來影響POWER。例如浅蚪，重復(fù)測量設(shè)計(jì)實(shí)際上總是比橫截面設(shè)計(jì)更強(qiáng)大藕帜；與較少的數(shù)據(jù)收集相比，響應(yīng)變量的更多時(shí)間點(diǎn)也會增加功率惜傲。效應(yīng)大小較小的兩組通常需要比具有更大效果大小的那些更大的功率來檢測洽故。統(tǒng)計(jì)程序和模型也會影響功效：當(dāng)模型的假設(shè)被違反時(shí)，非參數(shù)模型可能比參數(shù)模型更強(qiáng)大盗誊。相互作用項(xiàng)通常需要比主效應(yīng)更強(qiáng)的能量來檢測时甚。有時(shí)，時(shí)間點(diǎn)之間相關(guān)性也會影響功率哈踱。響應(yīng)變量的對數(shù)變換或?qū)B續(xù)響應(yīng)變量的歸類也會導(dǎo)致能量的損失荒适，因?yàn)檎{(diào)整模型會失去更多的能量，而Logistic或有序Logistic回歸往往比普通的最小二乘回歸需要更多的對象开镣。最后刀诬，通常情況下，丟失數(shù)據(jù)會降低功耗邪财，而不好的補(bǔ)償方法會極大地降低功耗陕壹。在某些統(tǒng)計(jì)程序中，個(gè)案刪除是處理丟失數(shù)據(jù)的默認(rèn)設(shè)置树埠。個(gè)案刪除是失去權(quán)力的最大因素之一糠馆。因此，當(dāng)您設(shè)計(jì)您的研究時(shí)弥奸，您需要盡一切可能將丟失的數(shù)據(jù)降至最低榨惠。
5.2 多樣性差異使用T檢驗(yàn)功效分析
①連續(xù)結(jié)果的冪公式:在微生物群研究中，樣本中單個(gè)類群的數(shù)量和多樣性可以用指數(shù)度量來概括盛霎，如Shannon多樣性指數(shù)(均勻度)和Chao1指數(shù)(豐富度)赠橙。基本目標(biāo)之一是比較不同群體的物種多樣性愤炸。調(diào)查人員提出的問題是：一個(gè)群體是否比另一個(gè)群體存在更多的多樣性期揪。或者规个，來自兩個(gè)數(shù)據(jù)集的alpha多樣性的分析結(jié)果彼此一致凤薛？在統(tǒng)計(jì)假設(shè)框架下姓建，可以提出問題。H0：組1中的分集與組2中的分集沒有區(qū)別缤苫；Ha：組1的多樣性與組2的不同速兔。多樣性度量可以形成假設(shè)檢驗(yàn)、冪和樣本量計(jì)算的基礎(chǔ)活玲。假設(shè)檢驗(yàn)將在第7-12“章節(jié)”中介紹涣狗。重點(diǎn)介紹功效和樣本量的計(jì)算。零假設(shè)是H0：u1=u2舒憾，其中u1和u2分別是組1和組2的真實(shí)分集平均值镀钓。另一種假設(shè)是Ha：u1不等于u2。零和可選的假設(shè)可以重寫為H0：u1-u2=0和Ha：u1-u2= δ不等于0镀迂。功效分析的目標(biāo)是以高概率檢測大小δ的差異丁溅。一般功率公式由下式給出
對于連續(xù)端點(diǎn)，檢驗(yàn)統(tǒng)計(jì)量為兩組的平均差δ探遵，然后
冪取決于α窟赏、δ、σ和n别凤。如果σ已知饰序，則在計(jì)算中使用正態(tài)分布，如果需要估計(jì)σ规哪，則使用非中心t(或表)求豫。假設(shè)σ已知，且n1=n2=n诉稍，那么我們可以使用數(shù)據(jù)來修改假設(shè)捌臊，如下所示：
如果出現(xiàn)以下情況破停，則拒絕H0分發(fā)
檢驗(yàn)均數(shù)差異的一般樣本量公式如下：
如果σ未知并且需要估計(jì)端铛，并且假設(shè)n1=n2=n练般，那么我們可以將假設(shè)修改如下：如果出現(xiàn)以下情況，則拒絕H0分發(fā)
如果每組的樣本量不平衡服爷，n1不等于n2杜恰，我們可以推導(dǎo)出樣本量公式如下
②ALS研究的多樣性數(shù)據(jù)：為了說明使用t檢驗(yàn)檢驗(yàn)多樣性差異的能力分析，我們比較了從肌萎縮側(cè)索硬化癥(ALS)研究中計(jì)算的G93A轉(zhuǎn)基因小鼠的Shannon多樣性仍源，該研究是關(guān)于丁酸治療對糞便微生物群的影響心褐。G93A小鼠攜帶人類肌萎縮側(cè)索硬化癥(ALS)引起的SOD1突變，這些突變概括了肌萎縮側(cè)索硬化癥(ALS)患者的神經(jīng)元和肌肉損傷笼踩。這些小鼠被廣泛用于研究肌萎縮側(cè)索硬化癥的病理機(jī)制和試驗(yàn)治療(張等人逗爹。2017年)。丁酸處理組9例嚎于，未處理組G93A突變體7例掘而。計(jì)算了16只小鼠的香農(nóng)多樣性挟冠。圖5.1中的分布顯示，丁酸處理組產(chǎn)生更高的多樣性袍睡，因?yàn)榕c不使用丁酸的對照組相比知染，該組的直方圖向右移動(更高的多樣性值)(注：垂直虛線標(biāo)記每組的平均多樣性)。在此數(shù)據(jù)中斑胜，丁酸處理組的香農(nóng)多樣性的平均值和標(biāo)準(zhǔn)差分別為2.504和0.170持舆，而不加丁酸的對照組的香農(nóng)多樣性平均值和標(biāo)準(zhǔn)差分別為2.205和0.209。對于連續(xù)變量伪窖，效應(yīng)大小正好是這兩個(gè)均值的差值(δ=2.504-2.205)。
圖5.1是通過以下R代碼生成的居兆。首先覆山，我們加載分類群多度表，并將原始的按樣本分類單元格式表轉(zhuǎn)置為按分類單元分類單元格式表
> ##load abundance table
> abund_table=read.csv("ALSG93AGenus.csv",row.names=1,check.names=FALSE)
> abund_table_t<-t(abund_table)
然后泥栖，我們從vegan包中調(diào)用多樣性函數(shù)來計(jì)算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)))
第三，我們根據(jù)示例信息創(chuàng)建變量“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
...
整個(gè)數(shù)據(jù)集包括1個(gè)月魏割、3個(gè)月、3.5個(gè)月和4個(gè)月的樣本數(shù)據(jù)钢颂。我們感興趣的是3到3.5個(gè)月期間治療和不治療的比較钞它。因此，我們將3個(gè)月到3.5個(gè)月的數(shù)據(jù)子集如下
> 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
最后殊鞭，我們根據(jù)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函數(shù)power.t.test()計(jì)算冪或樣本量：為了檢驗(yàn)香農(nóng)多樣性中沒有差別的零假設(shè)，可以使用t檢驗(yàn)或Wilcoxon秩和檢驗(yàn)操灿。我們把多樣性的假設(shè)檢驗(yàn)留到第7章锯仪。這里，我們重點(diǎn)說明如何使用R軟件計(jì)算冪或樣本量趾盐。在R中庶喜，可以使用Basic R中的函數(shù)power.t.test()和pwr包中的函數(shù)pwr.t.test()進(jìn)行功率分析。我們使用power.t.test救鲤。此函數(shù)的用法如下所示：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是每組樣本大小的數(shù)目，delta 是均值的真差蜒简，sd是標(biāo)準(zhǔn)差瘸羡，sig.level是顯著性水平(類型I錯(cuò)誤概率)，功率是測試的功率(1減去類型II錯(cuò)誤概率)搓茬，類型是t測試的類型犹赖，并且備選是單側(cè)或雙側(cè)測試队他。由于均值差的標(biāo)準(zhǔn)差未知，因此需要使用(5.5)中的公式進(jìn)行估計(jì)峻村。
> 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
現(xiàn)在可以通過調(diào)用power.t.test()函數(shù)獲得功率麸折。
> 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小鼠的大小樣本粘昨，隨機(jī)分配到丁酸處理或不處理對照垢啼，將提供83%的功率來拒絕兩組香農(nóng)多樣性沒有差異的零假設(shè)。如果樣本量增加到每組3個(gè)张肾，則功率將增加到99%以上芭析。我們可以生成功率和樣本大小圖，以可視化我們使用以下R代碼拒絕零假設(shè)所需的功率和樣本大小吞瞪。
> 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")
實(shí)際上馁启，冪和樣本量的計(jì)算可以在幾乎任何統(tǒng)計(jì)軟件包中完成(例如，在R中芍秆，我們使用power.t.test(n=2：10惯疙，delta=2.504-2.205，sd=0.05012)來生成圖5.2的數(shù)據(jù))妖啥。假設(shè)比較組中的一個(gè)是預(yù)先指定的微生物群組霉颠，則需要單樣本t檢驗(yàn)來估計(jì)功率或樣本大小。函數(shù)power.t.test()的默認(rèn)測試是兩個(gè)樣本荆虱。您可以指定type=“one.sample“進(jìn)行單樣本功率分析蒿偎，如下所示。
> power.t.test(n=2:10,delta=2.504-2.205,sd=0.05012, type = "one.sample")
5.3使用方差分析比較兩個(gè)以上組的分集的功率分析
①單因素方差分析的假設(shè)和冪理論：對于涉及兩個(gè)以上組的比較酥郭，可以使用單向方差分析(ANOVA)。例如愿吹，要比較三個(gè)或更多組中的多樣性不从，假設(shè)可以是：H0：三個(gè)組的平均多樣性相等；Ha：在這三個(gè)群體中犁跪，并不是所有的多樣性都是平等的椿息。差分析采用F檢驗(yàn)。要在完全隨機(jī)設(shè)計(jì)中計(jì)算一個(gè)處理或條件在k個(gè)水平上的全局F檢驗(yàn)的功效坷衍，我們首先需要理解和定義假設(shè)寝优。方差分析的基本思想是將多樣性的總體方差劃分為反映處理組或條件(因素水平)之間的差異和[由于測量誤差(殘差)]的處理組或條件內(nèi)的差異的分量。對于出現(xiàn)在i1枫耳；…乏矾；K水平的因子a，每個(gè)水平有j1；…钻心；J個(gè)觀測值凄硼，典型的單因素方差分析模型可以表示為
根據(jù)統(tǒng)計(jì)假設(shè)框架，假設(shè)被定義為捷沸，
其中摊沉，ui=組i的平均值，k=組數(shù)痒给，j=實(shí)驗(yàn)單位说墨。均值相等的F檢驗(yàn)采用單因素方差分析(ANOVA)假設(shè)數(shù)據(jù)是正態(tài)的，具有共同的組方差苍柏。也是N大于等于K尼斧，ni大于等于1，其中N是總樣本量试吁，ni是第i組的樣本量突颊。在零假設(shè)下，F(xiàn)統(tǒng)計(jì)量的分布服從中心F分布潘悼，而在另一種假設(shè)下，F(xiàn)統(tǒng)計(jì)量的分布服從帶有非中心性參數(shù)k的非中心F分布爬橡。因此治唤，當(dāng)零假設(shè)為真時(shí)，它遵循中心F分布糙申；當(dāng)它為假時(shí)宾添，它遵循非中心F分布。因此柜裸，冪可以定義為F統(tǒng)計(jì)量服從非中心分布的概率缕陕。?確切的功率是這樣給出的，
上述功率用兩個(gè)假設(shè)定義了兩條分布曲線疙挺。零假設(shè)下的F分布定義了中心F分布扛邑，而在備選假設(shè)下的分布定義具有非中心性參數(shù)λ的相同F(xiàn)。R使用公式(5.10)計(jì)算冪铐然。非中心性參數(shù)k是該公式中的關(guān)鍵參數(shù)蔬崩。在R中λ是如何定義的？對于平衡設(shè)計(jì)搀暑，設(shè)n=平衡組的樣本大小沥阳，則λ定義為
②使用R函數(shù)pwr.vova.test()計(jì)算冪或樣本量：為了說明如何使用pwr.vova.test()獲得不同樣本大小的能力，我們使用以下過程自点。首先桐罕，我們將四組1個(gè)月和4個(gè)月的數(shù)據(jù)分成兩組，分別為治療組和對照組。
> 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統(tǒng)計(jì)量溅潜。
> 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包中調(diào)用pwr.anova.test()函數(shù)來計(jì)算功率死宣。
> 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
上述結(jié)果表明伟恶，根據(jù)本次先導(dǎo)研究中檢測到的效應(yīng)大小，采用方差分析檢驗(yàn)毅该，每組需要52個(gè)樣本才能獲得80%的效果博秫。
5.4 比較組間的分類單元
①比較比例的假設(shè)和基本冪和樣本量公式：在微生物組數(shù)據(jù)集中，0表示樣品中沒有分類單元眶掌，1表示樣品中存在分類單元挡育。這樣的二進(jìn)制數(shù)據(jù)可以容易地轉(zhuǎn)換成比例朴爬。如果研究對特定分類群的豐富度感興趣召噩，那么可以使用卡方檢驗(yàn)來比較不同類群之間的比例母赵“汲埃跨組比較單個(gè)指定分類單元的假設(shè)為?H0：組1中指定分類單元的百分比與組2中的沒有區(qū)別；Ha：組1中指定分類單元的百分比與組2中的不同构韵。方程式(5.4)可以很容易地修改以測試分類單元的比例凶朗。當(dāng)我們測試分類單元比例時(shí)棚愤，我們正在評估對一個(gè)組(例如遇八，treatment)響應(yīng)的比例P1是否不同于對另一個(gè)組(例如， control）響應(yīng)的比例P2斯够。這等價(jià)于零假設(shè)H0：P1-P2=0與另一種假說Ha：P1-P2=δ≠0读规。
一般來說铃在，當(dāng)結(jié)果是二進(jìn)制時(shí)定铜，所需的樣本大小如下所示。
類似地畸陡，在連續(xù)情況下，可以推導(dǎo)出兩組比例的不平衡樣本量公式牲览。如果為n1≠n2跛蛋，則兩組樣本大小之間的比率為r=n1÷n2赊级。測試統(tǒng)計(jì)數(shù)據(jù)構(gòu)造如下：
②使用R函數(shù)power.pro.test()進(jìn)行功效分析：在我們的ALS研究中發(fā)現(xiàn)，口服2%丁酸鈉2.5個(gè)月與丁酸產(chǎn)生菌Butyribrio fisolvens在物種水平上的豐度顯著提高有關(guān)兑徘。我們有興趣比較這個(gè)分類單元在處理組和對照組之間是否以不同的速率存在。為了設(shè)計(jì)未來的研究崭闲，我們希望計(jì)算樣本量橄仍，以確保研究在此先導(dǎo)研究的基礎(chǔ)上有足夠的能力來檢測效應(yīng)大小侮繁。在本節(jié)中，我們將演示函數(shù)power.pro.test()斋射，以便使用我們的ALS G93A小鼠數(shù)據(jù)進(jìn)行功率分析，說明X2test和Fisher的精確測試腹躁。我們對丁狀弧菌的豐度計(jì)數(shù)數(shù)據(jù)進(jìn)行了轉(zhuǎn)換將纖維隔離物轉(zhuǎn)化為二元變量哑了，以指示是否存在纖維隔離物丁狀核弧菌。我們在表5.2中總結(jié)了轉(zhuǎn)換后的數(shù)據(jù)拆火。
下面的R代碼加載豐度表
abund_table_Spe=read.csv("ALSG93AButyrivibrioSpecies.csv",row.names=1,check.names=FALSE)
原始豐度數(shù)據(jù)是一個(gè)以物種為行、以樣本為列的計(jì)數(shù)表。下面的轉(zhuǎn)置函數(shù)“t()”將該物種計(jì)數(shù)表轉(zhuǎn)換成行的樣本和列的物種胞皱。
> abund_table_Spe<-t(abund_table_Spe)
可以從樣本標(biāo)識符獲得分組信息雾鬼。
> 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個(gè)樣本隨機(jī)分為丁酸治療組务甥，另外7個(gè)樣本設(shè)為對照組态辛。以下R代碼用于創(chuàng)建一個(gè)組變量來表示分組信息奏黑。
> grouping$Group <-with(grouping,ifelse(as.factor(X1)%in%c("B11","B12","B13","B14","B15","D11","D12","D13","D14"),c("Butyrate"), c("Control")))
將物種豐度表和分組信息表合并成一個(gè)數(shù)據(jù)框后馁害，即可進(jìn)行數(shù)據(jù)分析。
> 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的丁狀弧菌計(jì)數(shù)編碼為“不存在”，否則編碼為“存在”计雌。
> 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
數(shù)據(jù)匯總為2×2列聯(lián)表鼠渺。
> 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個(gè)對照樣品中有4個(gè)(57%)含有丁酸弧菌的比例鹃祖，而9個(gè)丁酸鹽樣品中有9個(gè)(100%)含有丁酸治療樣本做到了。使用(5.16)和(5.17)中的上述公式祖能，使用以下純R碼來實(shí)現(xiàn)樣本大小和功效計(jì)算。
> 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函數(shù)power.pro.test鳞滨。您可以指定多個(gè)樣本來測試功效澡匪。
>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
結(jié)果表明，在我們前期研究的基礎(chǔ)上荠瘪，每組11個(gè)樣本可以獲得83%的能量來檢測效應(yīng)大小。
③基于X2檢驗(yàn)和Fisher精確檢驗(yàn)的功率分析
比較比例的X2檢驗(yàn)的功率理論：檢驗(yàn)微生物群組成比例的假設(shè)也可以用卡方統(tǒng)計(jì)檢驗(yàn)。X2test統(tǒng)計(jì)量被定義為吠各，
與F統(tǒng)計(jì)量的分布類似藕筋，當(dāng)零假設(shè)為真時(shí)隐圾，X2統(tǒng)計(jì)量的分布服從中心卡方分布暇藏。當(dāng)它為假時(shí)把兔，它遵循具有非中心性參數(shù)λ的非中心卡方分布。因此某饰，本質(zhì)上诫尽，卡方分布的冪是數(shù)據(jù)來自非中心卡方分布的概率剂跟〔芮ⅲ基本上，X2分布定義了兩條曲線：一條是零假設(shè)下的中心X2分布，另一條是帶有非中心性參數(shù)k的非中心X2分布阐斜，零假設(shè)下概率為0.05的卡方用垂直線表示到推。如果X2統(tǒng)計(jì)數(shù)據(jù)落在行的左側(cè)惕澎，則表明它來自空分布捣卤，如果它在行的右側(cè)，則它來自替代分布子姜。線右側(cè)的備選假設(shè)曲線下的面積表示測試的威力牧抽。例如，在微生物組研究中，我們假設(shè)組1或條件1的比例遵循中心X2分布晨炕，而組2或條件2遵循非中心X2分布。
Using the Function pwr.chisq.test()實(shí)現(xiàn)功效分析：為了找到指定功率所需的樣本大小，R找到滿足公式的N。在這一部分中，我們使用我們的丁狀弧菌纖維隔離物存在/不存在的數(shù)據(jù)戈二，使用來自ALS G93A小鼠數(shù)據(jù)集的數(shù)據(jù)，說明使用v2test和Fisher精確測試進(jìn)行功率分析。要進(jìn)行功率分析，我們首先需要估計(jì)影響的大小嗓蘑。對于連續(xù)測度豌汇，例如在香農(nóng)分集的情況下梅尤，效應(yīng)大小僅計(jì)算為均值的差值。在分類數(shù)據(jù)中，有幾個(gè)度量可以用來作為效應(yīng)大小钝腺，包括Cramer‘s V(φ系數(shù)u)、優(yōu)勢比和相對風(fēng)險(xiǎn)(Cohen，1988)。你可以根據(jù)你的學(xué)習(xí)興趣選擇一個(gè)刮便。這里我們使用Cramer‘s V，它是名義變量關(guān)聯(lián)性的度量。實(shí)際上入客，它是重新調(diào)整的皮爾遜卡方統(tǒng)計(jì)數(shù)據(jù)，其值在0到1之間卓舵。
其中肿嘲，X2是零假設(shè)為真時(shí)的皮爾遜卡方，N是總樣本大小雳窟，t是行數(shù)減一或列數(shù)減一中的較小者尊浪。如果r是行數(shù)，c是列數(shù)封救，則t=最小值(r?1，c?1)誉结。當(dāng)然工育，對于2X2列聯(lián)表，這只是卡方除以觀測次數(shù)的平方根搓彻，等于φ系數(shù)u如绸。Cramer‘s V可以通過使用LSR軟件包中的函數(shù)CramersV來估計(jì)。
> library(lsr)
> cramersV(tbl)
[1] 0.3833
給定Cramer‘s V旭贬，現(xiàn)在我們可以使用PWR包中的函數(shù)pwr.chisq.test()進(jìn)行功率分析怔接。
> 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
結(jié)果表明，每組需要54個(gè)樣本才能在X2檢驗(yàn)的基礎(chǔ)上稀轨，以Cramer‘s V(φ系數(shù)u)為效應(yīng)大小扼脐，以80%的權(quán)值正確拒絕零假設(shè)。
Using the Functions
power.fisher.test() and power.exact.test()進(jìn)行功效分析：如果列聯(lián)表中的單元格值較小(<5)奋刽，則卡方檢驗(yàn)可能不正確瓦侮，則應(yīng)用Fisher精確檢驗(yàn)。因此佣谐，我們將分別從statmod和Exact包中說明最近開發(fā)的兩個(gè)R冪函數(shù)power.fisher.test()和power.exact.test()肚吏。與pwr.chisq.test()相比，兩個(gè)版本的Fisher測試在一次測試中指定多個(gè)樣本的靈活性較低狭魂。在下面的示例中罚攀，我們使用列聯(lián)表中的相同比例(p1=1.0党觅，p2=0.57)和每組15個(gè)樣本來估計(jì)功率。
> library(statmod)
> power.fisher.test(1.0,0.57,15,15,alpha=0.05, nsim=1000)
[1] 0.844
以下R代碼使用power.exact.test函數(shù)計(jì)算無條件精確測試的功率斋泄。
> library(Exact)
> power.exact.test(1.0, 0.57, 15, 15, method="Fisher")
$power
[1] 0.8454
$alternative
[1] "two.sided"
$method
[1] "fisher"
結(jié)果表明杯瞻，每組需要15個(gè)樣本才能在Fisher精確檢驗(yàn)的基礎(chǔ)上以84%的權(quán)值正確拒絕零假設(shè)。
5.5 利用Dirichlet-Polyomial模型比較所有類群的分布頻率
①多元假設(shè)檢驗(yàn)與Dirichlet-多項(xiàng)式模型：上面比較感興趣的分類單元的方法是單變量的“一次一個(gè)分類單元”分析炫掐。由于微生物組數(shù)據(jù)具有多變量結(jié)構(gòu)魁莉，我們經(jīng)常通過跨組比較每個(gè)類群的豐度來進(jìn)行多個(gè)類群的比較，然后進(jìn)行調(diào)整以進(jìn)行多次比較募胃。單變量方法通常不如多變量方法強(qiáng)大旗唁，因?yàn)樗鼪]有考慮到分類群之間存在的相互作用，而多變量方法聯(lián)合建模微生物類群豐度摔认。微生物組研究人員一直在努力開發(fā)統(tǒng)計(jì)方法逆皮，以進(jìn)行關(guān)于處理或?qū)嶒?yàn)因素對整個(gè)細(xì)菌分類群組合的影響的多變量假設(shè)宅粥，并估計(jì)此類實(shí)驗(yàn)的樣本量参袱。在比較微生物群的多變量假設(shè)檢驗(yàn)框架下，無效假設(shè)之一可以是：H0：不同的處理或?qū)嶒?yàn)因素對整個(gè)類群的影響沒有不同秽梅。Ha：不同的處理或?qū)嶒?yàn)因素對整個(gè)類群的影響是不同的抹蚀。多變量統(tǒng)計(jì)方法之一是Dirichlet多項(xiàng)式分布，它已被證明能很好地模擬微生物組數(shù)據(jù)企垦。La Rosa及其同事對Dirichlet多項(xiàng)式模型進(jìn)行了重新參數(shù)化环壤，使其適用于基于位置(均值比較)和尺度(方差比較/離散度)之間的差異進(jìn)行跨組假設(shè)檢驗(yàn)。該方法具有參數(shù)估計(jì)钞诡、多變量假設(shè)檢驗(yàn)郑现、冪和樣本量計(jì)算等功能。在這一章中荧降，我們首先介紹了重新參數(shù)化的Dirichlet多項(xiàng)式模型接箫，然后重點(diǎn)介紹了使用HMP軟件包計(jì)算冪和樣本量。
這種結(jié)構(gòu)的計(jì)數(shù)數(shù)據(jù)一般采用多項(xiàng)式分布進(jìn)行分析朵诫。然而辛友，它使用多項(xiàng)分布的適宜性假設(shè)每個(gè)類別的真實(shí)頻率在所有樣本中都是相同的(La Rosa等人。(2012年)剪返。此外废累，當(dāng)數(shù)據(jù)呈現(xiàn)過度離散性時(shí)，多項(xiàng)式模型可能會導(dǎo)致I類誤差增加脱盲。因此邑滨，使用多項(xiàng)式分布來模擬微生物組計(jì)數(shù)數(shù)據(jù)是不合適的，因?yàn)橛捎跇颖镜目勺冃郧矗⑸锝M數(shù)據(jù)中的每個(gè)分類單元在所有樣本中并不相同驼修。為了將微生物群分類單元計(jì)數(shù)數(shù)據(jù)與該數(shù)據(jù)結(jié)構(gòu)相匹配殿遂，重新參數(shù)化了Dirichlet-多項(xiàng)分布。其特征在于兩組參數(shù)：π和⊙：
其中乙各，⊙是超色散參數(shù)墨礁，或色散的量度。它測量樣本內(nèi)相對于多項(xiàng)分布的變異性過剩耳峦；⊙≥0表示過度分散恩静。與多項(xiàng)分布相比，Dirichlet多項(xiàng)分布的參數(shù)化使得它適合于通過比較π向量和比較⊙值來進(jìn)行基于位置差異的跨組假設(shè)檢驗(yàn)蹲坷。例如驶乾，對于兩組比較(即對照和病例)，零假設(shè)：π1=π2表示群落組成相等循签。有人建議使用Koehler和Wilson廣義Wald檢驗(yàn)來檢驗(yàn)零假設(shè)级乐。作為Dirichlet-Polyomial，重新參數(shù)化的Dirichlet-Polyomial模型既可用于分析分類單元組成數(shù)據(jù)县匠，也可用于分析秩豐度分布(RAD)數(shù)據(jù)风科。這兩種方法分別稱為“分類群組成數(shù)據(jù)分析”和“秩豐度分布數(shù)據(jù)分析。兩種方法有不同的側(cè)重點(diǎn)：一種是群落組成(那里有什么細(xì)菌)乞旦；另一種是群落結(jié)構(gòu)(如豐富度和多樣性)贼穆。
②Dirichlet-多項(xiàng)式模型下的冪和樣本量計(jì)算：Dirichlet多項(xiàng)式模型的冪和樣本量取決于概率模型參數(shù)、假設(shè)正在測試中兰粉，以及效果的大小故痊。模擬研究表明，效應(yīng)大小玖姑、超分散和樣本量影響功率愕秫。功率隨效應(yīng)尺寸的增大或樣本量的增大而增大，隨過色散的減小而增大焰络。在一些示例中戴甩，在保持效應(yīng)大小、超色散和樣本大小恒定的情況下舔琅，讀取次數(shù)也影響功率等恐，功率隨著讀取次數(shù)的增加而增加”蛤荆總之课蔬，統(tǒng)計(jì)能力受到每個(gè)處理的樣本數(shù)量(樣本大小或復(fù)制水平)和每個(gè)樣本的序列讀取數(shù)量(測序深度)的影響，以及過度分散郊尝。為簡單起見二跋，讓我們考慮一下在兩個(gè)實(shí)驗(yàn)組之間比較分類群頻率的微生物群樣本的問題，例如用丁酸處理的小鼠和不用丁酸處理的小鼠流昏。為了將分類群頻率與先前指定的微生物種群扎即、多于兩個(gè)組以及其他假設(shè)檢驗(yàn)進(jìn)行比較吞获。對于兩個(gè)組之間的π的比較，假設(shè)每個(gè)組的模型參數(shù)π和⊙都是已知的谚鄙，則檢驗(yàn)假設(shè)形成為：H0:?π1=π2 versus the alternative HA:?π1≠π2
影響的大小由分類頻率π1和π2的矢量彼此之間的距離來定義各拷。有幾種方法可以測量效果大小。其中之一是(5.20)中的Cramer‘s V闷营。對于表5.3中的微生物組數(shù)據(jù)結(jié)構(gòu)烤黍，行r是要比較的組的數(shù)量，列c是分類群的總數(shù)傻盟。根據(jù)公式速蕊，Cramer‘s V取決于樣本大小和讀取次數(shù)。為了減少樣本大小和讀取次數(shù)的影響娘赴，提出了一個(gè)修正的Cramer u準(zhǔn)則规哲，范圍從0到1，前者表示兩個(gè)類群的分類單元頻率相同诽表，后者表示分類單元頻率差異最大唉锌。
③使用HMP軟件包計(jì)算功率和尺寸：HMP是一個(gè)R包，用于假設(shè)檢驗(yàn)和能力計(jì)算关顷，用于比較最初來自HMP的元基因組樣本(La Rosa等人糊秆。2016)武福。功率和樣本量的計(jì)算可以通過使用統(tǒng)計(jì)量的漸近分布或通過HMP軟件包實(shí)現(xiàn)的Monte-Carlo模擬來進(jìn)行议双。在這里伴郁，我們使用ALS G93A小鼠的例子汽纤，通過蒙特卡羅模擬說明功率和樣本大小的計(jì)算狠鸳。
準(zhǔn)備數(shù)據(jù)集以使用HMP包：首先志衣，安裝并調(diào)用HMP包
> install.packages("HMP",repo="http://cran.r-project.org", dep=TRUE)
> library(HMP)
然后設(shè)置工作目錄耿戚，加載類群豐度數(shù)據(jù)集
> 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)
原始數(shù)據(jù)集具有按樣本分類的格式铃绒。如表5.3所示插爹，重新參數(shù)化的Dirichlet多項(xiàng)式模型和HMP包需要數(shù)據(jù)集具有分類抽樣格式别伏。因此莹规，我們使用t()函數(shù)將原始數(shù)據(jù)集轉(zhuǎn)置為具有Sample by Taxa格式的數(shù)據(jù)集赔蒲。
> 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
以下代碼檢查每個(gè)數(shù)據(jù)集中的分類群數(shù)和樣本數(shù)
> 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
利用分類群組成數(shù)據(jù)分析計(jì)算功率和大小：函數(shù)MC.Xdc.Statistics()使用模擬和似然比檢驗(yàn)來提供幾個(gè)樣本Dirichlet-多項(xiàng)式參數(shù)檢驗(yàn)比較的冪和大小良漱。這里舞虱，我們比較“Buty”和“NOButy”的兩個(gè)樣本集。一個(gè)示例語法如下：
MC.Xdc.statistics(group.Nrs, numMC = 1000, alphap, type = “ha”, siglev =0.05, est = “mom”)
其中母市，group.Nrs用于指定組中每個(gè)樣本的讀取次數(shù)/序列深度矾兜。NumMC用于指定蒙特卡羅實(shí)驗(yàn)的次數(shù)。在實(shí)踐中患久，至少應(yīng)該指定1000個(gè)椅寺。如果alphap用于計(jì)算測試統(tǒng)計(jì)數(shù)據(jù)的大小(類型I錯(cuò)誤)浑槽，即當(dāng)指定type=“hnull”時(shí)，則alphap指定矩陣返帕，其中行是參考組的alpha參數(shù)的向量桐玻。例如，以下代碼：Alphap<Fit_Buty$Gamma指定參考組的“Buty”伽馬矩陣荆萤。如果alphap用于測試統(tǒng)計(jì)數(shù)據(jù)的計(jì)算能力(類型II錯(cuò)誤)畸冲，即默認(rèn)值，或者當(dāng)指定type=“ha”時(shí)观腊，則alphap表示由每個(gè)組中每個(gè)分類單元的alpha參數(shù)矢量組成的矩陣邑闲。例如，以下代碼：alphap<-rbind(Fit_Buty$Gamma梧油，F(xiàn)it_NOButy$Gamma)為每個(gè)組中的每個(gè)分類群指定“Buty”和“NOButy”伽馬矩陣苫耸。如果type=“hnull”，則計(jì)算測試統(tǒng)計(jì)信息的大欣茉伞褪子；如果type=“ha”，則計(jì)算測試統(tǒng)計(jì)信息的冪骗村，這也是默認(rèn)設(shè)置嫌褪。Siglev用于指定測試或功率計(jì)算大小的重要級別。默認(rèn)值為0.05胚股。est用于指定要與似然比檢驗(yàn)統(tǒng)計(jì)量“mle”或“mom”一起使用的參數(shù)估計(jì)器的類型笼痛。默認(rèn)值為“mom”。HMP軟件包的作者指出琅拌，“最大似然估計(jì)”的運(yùn)行時(shí)間要長得多缨伊，并且在小樣本情況下不是最優(yōu)的，而“mom”的結(jié)果在小樣本情況下更為保守进宝。首先刻坊，使用函數(shù)DM.MoM()獲取數(shù)據(jù)的Dirichlet多項(xiàng)式參數(shù)列表(即loglik、Gamma党晋、pi和theta)谭胚。
> fit_Buty <- DM.MoM(Buty_t);fit_NOButy <- DM.MoM(NOButy_t)
> fit_Buty
其次，設(shè)置蒙特卡羅實(shí)驗(yàn)的次數(shù)未玻，這里我們使用1000灾而，這是實(shí)踐中推薦的最小值，并生成每個(gè)樣本的讀取次數(shù)深胳。
> 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)
第三绰疤，計(jì)算測試統(tǒng)計(jì)數(shù)據(jù)的大小(I類錯(cuò)誤)。
> alphap <- fit_Buty $gamma
> pval1 <- MC.Xdc.statistics(group_Nrs, numMC, alphap, ”hnull”)
> pval1
[1] 0.000999
最后舞终，測試統(tǒng)計(jì)數(shù)據(jù)的計(jì)算能力(類型II錯(cuò)誤)轻庆。
> 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樣本數(shù)據(jù)集中獲得的Dirichlet多項(xiàng)式參數(shù)癣猾，使用5%顯著性水平比較Buty和NOButy種群分類群頻率的功率分析。每個(gè)條目表示指定采樣數(shù)和0.05顯著性水平下的讀取數(shù)所獲得的功率余爆。例如纷宇，對于樣本數(shù)=10，并且每個(gè)樣本的讀取次數(shù)=500蛾方，研究具有98.8%的能力來檢測在數(shù)據(jù)中觀察到的效應(yīng)大小像捶。結(jié)果顯示，當(dāng)樣本較小時(shí)桩砰，增加讀取次數(shù)(在本例中為5)會影響功率拓春；但當(dāng)達(dá)到較大樣本大小時(shí)(在此情況下，為10)亚隅，增加讀取次數(shù)不會影響功率硼莽。結(jié)果還表明，當(dāng)獲得足夠的樣本時(shí)煮纵，增加樣本量并不會增加功率懂鸵。
利用秩豐度分布數(shù)據(jù)分析計(jì)算功率和尺寸：幾種樣本RAD-概率均值檢驗(yàn)與已知比例參考向量的比較
函數(shù)MC.Xmc.Statistics()使用模擬和廣義Wald-type統(tǒng)計(jì)量來提供幾個(gè)樣本RAD概率均值測試與已知比例參考向量比較的能力和大小。這里行疏，我們再次比較“Buty”和“NOButy”的兩個(gè)樣例集匆光。一個(gè)示例語法如下所示：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()函數(shù)中定義货葬。Pi0是RAD概率平均向量采幌。如果group.pi=“hnull”劲够，則忽略此參數(shù)震桶；如果group.pi=“ha”，則指定一個(gè)矩陣征绎，其中每行都是每個(gè)組的向量pi值蹲姐。Group.theta是每個(gè)組的超離散值的向量。分析稀有類群的技術(shù)挑戰(zhàn)是匯聚率低人柿，測試統(tǒng)計(jì)不精確柴墩。為了提高收斂速度和準(zhǔn)確估計(jì)，在微生物群研究中凫岖，一種方法是去掉頻率較低的分類群江咳，例如在任何百分比小于1的樣本中；另一種方法是將所有頻率較低的分類群合并起來哥放。分類單元歼指，例如爹土，將低于1%的所有類群加權(quán)平均為一個(gè)分類單元，稱為“其他”分類單元踩身。在上述分類群組成數(shù)據(jù)分析中胀茵，我們在調(diào)用MC.Xdc.Statistics()函數(shù)時(shí)沒有使用這些技術(shù)。在這里挟阻，我們使用Data.filter()函數(shù)按豐度遞減的順序?qū)Ψ诸悊卧M(jìn)行排序琼娘，并將不太豐富的分類單元折疊成一個(gè)標(biāo)記為“Other”的類別。一個(gè)示例語法如下所示：Data.filter(data, order.type = “sample”, minReads = 1000, numTaxa = 10,perTaxa = NULL)附鸽。其中脱拼，data是每個(gè)樣本(行)的分類計(jì)數(shù)(列)矩陣。如果order.type=“sample”坷备，則根據(jù)分類頻率對分類進(jìn)行排名挪拟；如果order.type=“data”，則根據(jù)所有樣本的累計(jì)分類計(jì)數(shù)對分類進(jìn)行排名(默認(rèn))击你。假設(shè)minReads=一個(gè)讀取截止值玉组，則讀取總數(shù)小于該讀取截止值的樣本將被刪除。參數(shù)numTaxa和perTaxa丁侄，一次調(diào)用中只能指定一個(gè)惯雳。參數(shù)numTaxa用于指定要保留的分類群數(shù)量，同時(shí)將其他(不太豐富的)分類群折疊為“其他”類別鸿摇。參數(shù)perTaxa用于合并要保留的數(shù)據(jù)的百分比石景，同時(shí)折疊剩余的分類群。首先拙吉，使用Data.filter()函數(shù)按豐度遞減的順序?qū)Ψ诸惾哼M(jìn)行排序潮孽，并將豐度較低的分類群折疊為“其他”類別。
> 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
其次筷黔，使用函數(shù)DM.MoM()獲取數(shù)據(jù)的Dirichlet多項(xiàng)式參數(shù)列表(即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
第三椎例，設(shè)置蒙特卡羅實(shí)驗(yàn)的次數(shù)，這里我們使用1000请祖，這是實(shí)踐中推薦的最小值订歪，并生成每個(gè)樣本的讀取次數(shù)。
> 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)
第四肆捕，設(shè)置各類群的類群頻率向量(類群比例)和超分散參數(shù)的值刷晋。
> pi0 <- fit_Buty$pi
> group_theta <- c(0.007523, 0.01615)
第五，計(jì)算測試統(tǒng)計(jì)數(shù)據(jù)的大小(I類錯(cuò)誤)。
> pval1 <- MC.Xmc.statistics(group_Nrs, numMC, pi0, group.theta=group_theta, type="hnull")
> pval1
[1] 0.08492
最后眼虱，計(jì)算測試統(tǒng)計(jì)數(shù)據(jù)的功率(類型II錯(cuò)誤)或舞。
> 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顯示了使用函數(shù)MC.Xmc.Statistics()進(jìn)行秩豐富分布(RAD)數(shù)據(jù)分析的功率分析結(jié)果。與分類群組成數(shù)據(jù)分析相比蒙幻，RAD數(shù)據(jù)分析獲得了更高的檢測效果大小的能力映凳。例如，每組僅5個(gè)樣本邮破，每個(gè)樣本500個(gè)讀取诈豌，我們可以達(dá)到89.61%的功率。實(shí)際上抒和，與分類群組成數(shù)據(jù)分析相比矫渔，RAD數(shù)據(jù)分析方法降低了功耗，因?yàn)楫?dāng)忽略分類群標(biāo)簽時(shí)摧莽，數(shù)據(jù)中的信息會丟失庙洼。在這種情況下獲得的更高的能力是因?yàn)镽AD數(shù)據(jù)分析將196個(gè)分類群減少到11個(gè)分類群，其中包括1個(gè)“不那么頻繁地匯集在一起的豐富的分類群”镊辕。
幾種樣本比例向量未知的RAD-概率均值檢驗(yàn)比較：函數(shù)MC.Xmcupo.Statistics()使用模擬和廣義的Wald-type統(tǒng)計(jì)量來提供幾個(gè)樣本RAD概率均值測試比較的功率和大小油够，而不需要比例參考向量。這里征懈，我們再次比較“Buty”和“NOButy”的兩個(gè)樣例集石咬。一個(gè)示例語法如下所示：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在函數(shù)MC.Xmc.Statistics()中定義维贺。函數(shù)MC.Xmcupo.Statistics()的運(yùn)行方式與函數(shù)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顯示了使用函數(shù)MC.Xmcupo.Statistics()進(jìn)行秩豐富分布(RAD)數(shù)據(jù)分析的功率分析結(jié)果幸缕。函數(shù)MC.Xmcupo.Statistics()實(shí)現(xiàn)的冪小于函數(shù)MC.Xmc.Statistics()的冪群发，這是因?yàn)闆]有使用比例的參考向量。
④使用HMP程序包計(jì)算效果：HMP包具有計(jì)算效果大小的功能发乔。這里，我們使用函數(shù)Xmcupo.ffectsize()來說明效果大小的計(jì)算雪猪，該函數(shù)計(jì)算測試統(tǒng)計(jì)量Xmcupo.sevsample()的Cramer‘s Phi和修改后的Cramer’s Phi Criteria栏尚。使用函數(shù)Xmcupo.sevsample()進(jìn)行假設(shè)測試。在這里只恨，我們將重點(diǎn)放在效果大小的計(jì)算上译仗。語法為：Xmcupo.effectsize(group.data)其中抬虽，group.data是一個(gè)列表，其中每個(gè)元素都是每個(gè)樣本(行)的分類計(jì)數(shù)(列)矩陣纵菌。
> ##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
在基于檢驗(yàn)統(tǒng)計(jì)量Xmcupo.sev樣本的RAD數(shù)據(jù)分析中觀察到的效應(yīng)大小分別為：Cramerφ0.03阐污，Modified-Cramerφ0.15。
5.6?總結(jié)
本章重點(diǎn)介紹微生物組數(shù)據(jù)的功率和樣本量計(jì)算咱圆。介紹了單變量和多變量分析方法笛辟。對于單變量分析方法，我們介紹并說明了參數(shù)和非參數(shù)檢驗(yàn)和模型序苏；對于多變量分析方法手幢，重點(diǎn)介紹了使用參數(shù)Dirichlet多項(xiàng)式模型計(jì)算冪和樣本量。在每個(gè)分析中都引入了統(tǒng)計(jì)理論和相關(guān)的R軟件包忱详。我們從介紹假設(shè)檢驗(yàn)和力量分析的要素開始围来，并利用微生物組數(shù)據(jù)形成假設(shè)檢驗(yàn)和力量分析。其次匈睁，我們劃分了四個(gè)部分监透，分別涵蓋力量分析和樣本量分析，分別使用t-test檢驗(yàn)多樣性差異航唆；使用方差分析比較兩個(gè)以上組的多樣性才漆；比較組間感興趣的分類群；以及使用Dirichlet-Polyomial模型比較組間所有類群的頻率佛点。微生物組能力分析對于微生物組研究的研究設(shè)計(jì)是非常重要的醇滥。文獻(xiàn)中的假設(shè)檢驗(yàn)方法和模型有限，特別是針對合適的多變量模型的假設(shè)檢驗(yàn)更是鳳毛麟角超营。
