统计04回：族群与采样

当我们碰到一大群数字时，或是你想研究自然界的状况，你都不可能一次把所有的数据都量到手。想想看如果你想知道学校操场的草长得多好每棵多重，你一定不会把每一棵草都拿来称重。全校的草就是你希望研究的族群，而你唯一可行的做法就是找一个具有代表性的区块，把这区的草拿出来称重。这时，这一区的草就是你的样本。

基本上的想法是，我们打算用样本里的数字特性来推算估计族群里的数字特性。我们被迫要这样做，因为没办法把整个族群拿来研究。接下来当然就是希望知道样本里的数字和族群里的数字在特性上有没有什么不同？前面提过我们关心的是集中在哪里和分散程度，换成术语也就是平均值和标准差。样本里的平均值和标准差，和族群的平均值和标准差，有没有什么不一样？

先假设我们的采样没有问题，不管是随机采样还是系统采样，反正就是采好了。

我们希望知道族群平均值与样本的平均值会不会一样。这一点你用常理去想就知道会是一样的。你随便拿起来一颗苹果，用这课的重量猜一整箱的平均值，误差会很大，如果拿三颗的平均值来猜，受到特别大特别小的苹果的影响会被正常的苹果抵消，所以准一点。如果你拿半箱的平均值来猜，就更准了。如果你拿第一组（假设3颗好了）猜，发现低估了，拿两组，发现高估了，这样算下来，有的高估了，有的低估了，终究整箱猜完后高估低估会抵消。这不就是在算平均值吗？所以样本平均值的平均值等于族群的平均值，我们拿样本平均值来猜族群平均值的想法应该是合理的。可是你每拿三颗跟每组拿6颗的结果不太一样；每组6颗的状况下，各组平均值会比较接近族群平均值。想象一下有600位恐怖不受控制充满能量的幼稚园小朋友；每20个人一件教室坐好，跟大家一起坐草地上完全没人管的状况下，小朋友的能量是一样的。分教室后每一间都有很吵的跟安静的人，但是平均吵醒程度在各个教师差不多。我们用每间教师里的分贝数来推估全校平均吵闹程度，应该是合理的。但是这间教师里最吵的通常不是全最吵的，最安静的也不是全校最安静的人，拿全校来看，个体间的变异更大，一件教室只有一部分的变异存在。

于是你发现了，每个样本的平均值理论上会跟整个族群的平均值一样。而样本里的变异，只有族群里变异的一部分。

统计05回：族群分布。为什么要管分布？

过去数学界先贤先烈花了好大工夫找出公式来描述数字分布，如果有一群数字符合什么状况的时候，它们就符合什么样的分布，而且可以用公式来把它们找出来。这一套东西，和我们有什么关系？

当我们要研究自然界的状况时，不可能把所有的数值都收集到，所以我们要采样。我发现拿到的麦当劳薯条好像有点轻，我怎么知道加薯条不容易控制数量本来就有多有少，还是我这包缩水了？如果我可以买10000包正常品规格的薯条，每一包都去测量重量，我就知道重量有多又少，特别多或特别少机会都不打。这就是薯条重量分布的状况，利用它，我们可以知道如果从这10000包里拿出任何一包薯条，重量就会比平均值多两克的机会有多少。

接下来我想知道手上这包疑似缩水的薯条到底是不是缩水了。有了分布图，那我只要比对一下，看看手上这包的重量出现的几率有多少。如果机会高，那就是正常误差范围，

不幸的是，通常我们不会去买10000包来证实手上这包的问题。所以我们只好从各种证据里去推论我们想猜的族群，它的平均值到底是用什么方法分布的。

统计06回：常用的分布

最常用的常态分布（normal distribution）。因很多分布的平均值都是常态分布，所以（打算盘）如果我们可以由很多样本分别计算平均值，这些平均值都是常态分布，所以就可以拿常态分布的公式来推算我们猜对猜错的几率了啊。

另一在生物上会用的分布是二项式分布（binomial distribution）。这个分布描述的是自然界中只有两种状况的数值分布，是与非，生与死等。丢铜板最容易理解了。你丢一个铜板出现正面的机会有多少？丢两个呢？三个？一百个？描述铜板个数和出现正面个数的关系，就要用到二项式分布。

统计07回：常态分布的标准化

怎么把一堆数字标准化？因为每个族群的数字集中在哪里都不一样，所以我们硬性规定标准化的第一步是大家都往0集中。但是分散程度还是不一样，所以也要针对标准差调整。

统计08回：样本平均值与族群数字的特性如何互换

前面提過, 為了要能用樣本計算出來的平均值推估族群的平均值, 我們必須要知道 樣本平均值 的分佈狀況. 現在起請把"樣本平均數"當做一個新的參數, 來討論一下這個新的參數到底有什麼特性

我們依然關心平均值和分散狀況. 假設樣本平均值的分佈符合常態分佈. 前面提到過, 不管你的數值是符合哪種分佈, 在大部份的狀況下, 平均值都是符合常態分佈的. 那我們可以來看我們關心的 “樣本平均數" 的平均值及分散狀況.

那樣本平均數的平均值會落在哪裡? 答案是族群平均值. 所以我們才能用樣本平均數的平均值來推估族群平均值.

那样本平均数的分散状况呢？如果我们拿到一堆从每个样本计算得到的样本平均数，这些数字之间的变异一定比原来族群小，因为有一些变异被同一样本里的其他数字抵消掉了。我们把每个数字变成平均值加与平均值之间的差异。

统计09回：student's t distribution

前面提過, 如果我們從未知族群裡採樣, 這個樣本的平均值大致會符合 normal distribution. 這個想法沒錯, 但是有個沒講的前提是你的樣本必須是大的. 可是我們不是每次都能採個100個數值的樣本啊! 如果我們的樣本裡的數值變少, 那用它來來估計族群平均值就比較可能猜錯, 於是前人幫我們找出了當樣本小的時候平均值該有的分佈應是多少, 然後把它叫做 Student’s t distribution.

由於 t distribution 跟真實狀況比較接近, 所以利用 t distribution 來處理數據問題就比較合理了.

Normal distribution 配上個 Z 值來標準化, 那 t distribution 也有個對應的 t 值.

统计10回：t test怎么用

在使用 t test 之前, 我們先來搞清楚狀況.

你手上有一把數值, 而且你可以用這一把數值算出平均值及標準差. 可是這一把數值只是眾多樣本之一, 你的目的不是樣本, 而是要用在樣本上看到的特性來回推真正族群裡的特性.

我們的目標是要知道我們在樣本裡看到的數字是否符合一個已知的數字. 例如你線上買了 10 包大包的楊記地瓜酥, 上面的標示是每包 350 克, 你想知道老闆到底有沒有裝滿 350 克因為每次都覺得還沒滿足就見底了. 你有個目標數字 350, 裝地瓜酥一定有誤差不會每一包都剛剛好 350 克. 而你不知道平常生產時的誤差有多大. 你這 10 包平均起來重 355 克, 會不會是在誤差範圍裡? 所以你能做的是用你的 10 包間的誤差來估算生產線的誤差狀況(分佈), 然後看看在這樣的分佈下出現 355 克的機會大不大. 機會大, 那你拿到的是正常產品, 機會很小, 那你可以考慮竊笑賺到或是打電話跟老闆說你的包裝標示不實(咦).

因為樣本不大, 所以你找上了 t distribution 來描述樣本平均值的動態. 回想上次的公式:

统计14回：ANOVA各组间是否有差别

前面谈的是一个实验组和一个对照组间的比较，利用比较平均数或者比较变异来知道实验组的试验处理有没有造成影响。

接下来看如果有很多不同的处理要比较的时候该怎么办。例如，我们想知道使用iphone，HEC，Samsung, Motorola和不用手机的人，在每天使用脸书的时间长短上有没有差别。这个分析必须同时比较很多组，不再是一对一的比较了。

这时我们要动用的是方差分析，意思是把这些数字间的变异拿出来看看，看再各组内的变异（同厂牌手机）比较大，还是在不同组间（不同厂牌手机）的变异比较大。

第一步是要把变异区分成组内变异或组间变异。组内变异是同一组的数值间计算变异，每个数值减去平均值，然后平方，再相加，就是这一组的变异（变方）。组间变异也很好算，就是这一组和另一组的差异再平方，把这一组的平方和全体总平均相减得到这组平均和总平均的差别，这个数值平方后乘上这一组的n值（用平均算好后，当然要乘上组员人数才能代表整组变异）。每一组组平均和总平均差异的平方总和，就是组间的变异。你可以想象，当组间变异比较大，组内变异比较小的时候，各组应该是组组分开的，就像喜宴开始大家还一桌桌有秩序的做好，你看得出各组的区别。如果组内变异大，组间变异小，那就是各桌客人跑到别桌串门子小朋友冲来冲去的大混乱场面，分桌已经没有意义了。回到统计上来看，如果各组数字间大混乱，那这些组之间是分不开的，我只好说各组处理间没有显著差异。如果各组间有明显分开的趋势，我就该说各组间处理有显著差异。

可是要用什么标准来说组间差异大不大？这个简单，我们只要拿平均的组间变异去跟平均的组内变异比较就好了。计算上组间变异容易算，所有数值的总变异也容易算，所以把总变异算出来，减去组间变异，就是组内变异了。要怎么平均呢？有5组，自由度是5-1=4，平均变异是把组间变异除以4；总变异的自由度是总共数值数减1，把总自由度减组自由度，剩下来的自由度就给组内变异去分了。真正计算时我们也不用去担心这些该怎么算，app会搞懂。你只要会看组间变异和组内变异的比值算出来的F值出现的几率有多大就可以了。

未來這一類的分析還有很多地方要用到, 請開始這樣看這些數據. 我們把變異依據可能造成變異的原因一一做分析, 像在這裡我們認為不同組的處理可能是造成變異的原因, 所以把因為分組造成的變異拿出來, 跟不能被解釋的變異(在這裡是組內變異, 採樣誤差/個體差異造成的)來比較.未來在不同的驗設計裡我們可能會看到不只一種可能造成影響的因子陸續出現, 想法也是把因為各因子造成的變異都分離出來看, 用總變異減掉這些可以被解釋的變異後剩下的是不能被解釋的變異(error). 把每一種可以被解釋的變異除以不能被解釋的變異, 看看它們跟不能被解釋的變異相比是不是夠大. 如果夠大, 就可以認定這個因子有造成影響了.

主坐标分析Principal Coordinates analysis (PCoA)

這是微生物群聚分析裡經常出現的一種分析方法. 背後的想法是這樣的: 我们的目标是找出来自不同样本里的群聚组成，然后想要知道不同群聚在组成上像不像。第一步是用你喜欢/觉得合适的方法定义任两个群聚组成间的距离来代项它们之间有多像。接着把这些距离转成欧几里得距离，再用主成分分析PCA来降维。所以第一步是找出距离，然后交给PCoA用这些距离重新转换成我们能理解的距离来呈现。

为什么有那么多种距离？想象你要去成大玩，想知道从慈大到成大有多远。如果我用google map量直线距离给你，你不会满意的，因为你不可能在中央山脉打个洞走直线过去。画一条直线的距离是欧几里得距离，但是我们有兴趣的不见得是这个距离。由于PCoA处理的是距离矩阵，不管你用的是什么古怪方法算出来的距离，都可以喂给PCoA,但是PCA就得给它符合欧几里得距离的资料才能处理了。

PCoA拿到距离矩阵后，是先把这些数据重新整理，以欧几里得距离呈现，而且不限要动用到多少个变因，反正把所有资讯都转换出来就好了。理论上，如果你有n个群聚，用n-1个变因就一定能完整呈现所有资讯。所以PCoA用了大量的变因来把不管是怎么得来的距离，转换放进由n-1个轴形成的欧几里得空间里。

不過這成了另一個問題. 這麼多個軸的空間, 是遠遠超過我們大腦可以想像的狀態, 我們要如何從裡面看趨勢哩? 這時我們再找老朋友 PCA 幫忙, 把 n-1 度空間的資料降維作為二度或三度空間, 就可以拿來看趨勢了.

统计15回：微生物生态的数据

微生物生態的分析需要利用到多變值統計分析方法。先來談談會碰到的數據。

我們在微生物生態上碰到的數據跟一般生態學裡碰到的群聚數據一樣，大概會有兩類。你到每個地方採樣，回來分析後每個樣本會有一組菌相組成的結果。把每個樣點的結果放在一起，你的數據大概會長這樣: 每一行是一個樣本，每一列是一個 OTU，組成一個大表。我們可以把這種各樣點菌相組成的資料叫做 sites x species table。

除了知道菌相組成，你還會想知道是什麼原因造成了你看到的菌相組成。影響生物組成最有可能的原因是環境條件，所以你一定會想在那個樣點同時測量環境因子，例如土壤含水量、氣溫、硝酸根含量等等。所以每個樣本會有一組環境因子測量的結果。把每個樣點的結果放在一起，你的數據大概會長這樣: 每一行是一個樣本，每一列是一個環境因子，組成一個大表。我們可以把這種各樣點環境因子的資料叫做 sites x environmental parameters table。

接下來，我們可以分別對這兩種矩陣做分析，也可以把兩個合起來用。但是不管用哪一種，我們都得同時處理很多種菌株的數據或很多個環境變因的數據，這時，就是多變值縭計方法上場的時候了。我們會跟著 GUSTA ME 來談談微生物生態裡用的多變值統計方法。

统计16回：微生物生态数据标准化

你的數據矩陣裡有很多樣本，很多環境因子，這些可是得先經過處理才能開始做分析的。第一步，一定是檢查有沒有抄錯記錄錯而出現很怪的數字，例如大家的數值平均在 100 左右，有個記錄值卻是 1，這種數字就要回實驗記錄本裡看看原始資料，看看是不是在輸入電腦時打錯了。這一步很簡單，但是一定要檢查，不然後面分析可能會因為一個小錯誤而必須重新來過。

接下來是標準化。拿菌相組成來當例子，你拿到的原始數據是 NGS 分析後的序列數，A 樣本可能有 100000 條序列，數量最多的 OTU01 可能有 10000 條序列。B 樣本可能有 200000 條序列，裡面的 OTU01 可能有 15000 條序列。我可以說在 B 樣本裡 OTU01 的數量比較多嗎?

你可以說在 B 樣本裡我們讀到的 OTU01 序列的數量比較多，但是這個結論沒有意義，因為我們在意的是相對數量。在 A 樣本裡 OTU01 佔了 10000/100000 = 10%，但在 B 樣本裡 OTU01 佔了 15000/200000 = 7.5%，所以 OTU01 在 B 樣本裡的比例是比較低的。從這個例子裡你可以看出來，對每個樣本的序列總數做標準化，是必須要做的步驟，不然定序時得到序列數比較多的樣本就會造成比較大的影響了。

另一個需要標準化的地方是環境因子。你測量土壤 pH 時，得到的數值大概在 6.0-8.0 之間，差距大約在零點幾而已。如果是記錄氣溫，數值會在 10-30 之間，差距大約在個位數。如果記錄海拔高度，你的數值可能從幾百到上千，差距可能有好幾百。放進統計方法裡計算時，海拔高度數字大，一定比 pH 容易造成影響。所以這時就要對每個環境因子做個標準化。當用的標準化方法是生統學過的 Z score，每個數值先減掉這個因子的平均值，再除以標準差。

统计17回：简介GUSTA ME里提到的各类统计方法

多變值統計在微生物生態上會用到的方法大概有下面幾大類。

第一类方法的目的是用sites x species (或sites x environmental parameters)的数据，来看各样点间有多像，然后呈现出谁跟谁比较像，谁跟谁比较不像的结果。这类方法在GUSTA ME里是放在（dis）similarlity-based methods里，因为在计算上是要依靠算各样点间的相似或相异程度来达成的。

第二类方法是同时使用sites x species 和sites x environmental parameters 的数据，来看各样点在species组成间的相似能不能以environment来解释。概念上来说，就是在有sites x species- 和sites x environmental parameters数据的情况下，第一类方法一次只看其中一组，第二类方法则是同时看两组，并用其中一组来解释另一组。这类方法因为是强迫环境资料来解释菌相组成的相似程度，在GUSTA ME里是放在Constrained analyses里。

第三类方法是要找出区分样本的方法。例如你有从森林来的样本，有平地样本，也有都市里的样本。你想知道来自这三地的样本在菌相组成上差在哪里。这时你的目的不是看它们是否相似，而是要找出可以把它们区分开的方法。这类方法在GUSTA ME里是放在discrimination里，就跟字面上意思一样，是一群目的是要做区分的分析方法。

第四类方法是检验假说是要用的。大部分前面提到的方法都是让我们看出趋势，AB样点比较像，AC比较不像，但是这些结果都很主观。如果你想真正知道某些点之间的差异到底明不明显，你需要利用统计检验。这类方法在GUSTA ME里是放在Hypothesis testing里。

第五类方法使用另一种策略来看环境因子和菌株组成间的关系。因为环境因子跟菌相组成是一个大表，我们可以先用统计方法把环境因子缩成少少的变因，当做代表环境状态的指标。接着再用统计方法把菌相组成也浓缩成少少的生物指标，然后看生物指标的变化是不是会跟环境因子指标有关系。源自这样的想法的统计方法在GUSTA ME里是放在indirect gradient analysis 里面。

统计18回：(dis) similarity methods 简介

在微生物生态研究中，我们好奇是采了这么多个样本，到底哪些样本之间会比较像，哪些样本之间会不同。如果我们可以用个标准方法来计算两两样本在菌相间的相似程度，这个目的就可以达成。

可是，要用什么方法来计算相似程度呢？这时就要用特定的相似性指数了。可用的相似性指数很多，可以根据目的来选用。在微生物生态里常用来描述菌相组成相似度的指数包括Bray-Curtis similarity 或UniFrac distance

有了距离后就可以根据距离来把所有样本摆好找趋势了。我们有两种方法可以来看趋势。一种是ordination,把所有点放在同一个空间，这通常会是人脑容易理解的二度或三度空间，然后自己看谁跟谁比较像。另一类方法是clustering，是把样本根据相似度用个树形图来表示。

由于有了哪些样点比较接近的趋势，你一定会想知道这些看起来好像很像的菌相组成，到底是真的一样还是假的一样。这时就需要统计方法来告诉我们这些看起来很像的菌相组成，它们真的一样的机会有多少。我们可以根据样本间的相似度来计算机率，进行假设检验（hypothesis thesting）。

统计19回：Cluster analysis

Cluster analysis 可以根据分群的方法再区分为hierarchical cluster analysis 和non hierarchical cluster analysis （层序聚类分析和非分层聚类分析）

在Hierarchical clustering 的结果会是把所有样本分成两群，每一群再往下再分两群，知道所有样本都分开来，是个阶层性的系统。而分层级聚类则不会把所有的样本都以二分罚分开，而是让你指定要分多少群，再把样本们分成你要的那么多群。

Hierarchical cluster analysis

目前最常用的方法是 agglomerative clustering (以逐次聚合的方式)，操作時是先把最近的樣本抓在一群變一組，然後找剩下來各樣本裡誰跟這一組最像，把它拉進來變成這一組的成員。這個更新過的組可以再在剩下來各樣本裡找出跟這一組最像的樣本後拉進來，用這個方法反覆到所有的樣本都拉進來為止。

在這個方法裡，要先決定計算相似度的方法。第二件要做的，是決定怎麼計算一點和一組之間的距離。例如，A組裡有 A1, A2, A3，我想算 B 和 A 組間的距離。如果用 complete linkage 計算法，A 和 B 的距離是 A1-B、A2-B、A3-B 三段距離的平均值。如果是 single linkage 計算法，則是 A1、A2、A3 之間最接近 B 的那點跟 B 的距離。不同計算方式，當然可能會出現不同結果。

你可以參考 GUSTA ME 裡關於 Hierarchical clustering analysis 的說明。

Non-hierarchical cluster analysis

常用的方法是 K-mean clustering。你先指明要分多少群(K群)，軟體先指定 K 個中心點，依離各中心點的距離把每個點分成 K 組。接着根据分到各组的样点计算出新的中心点，然后再用新的中心点来分群。一直持续到结果稳定为止。

你可以參考 GUSTA ME 裡關於 Non-hierarchical clustering analysis 的說明。

统计20回：欧几里得距离

你熟悉的世界就是歐幾里德空間，在這裡畢氏定理是存在的。座標(Xa, Ya)的點 A 和原點的距離是 Xa 平方加 Ya 平方再開個根號，座標(Xb, Yb)的點 B 和原點的距離是 Xb 平方加 Yb 平方再開個根號，而點 A 和點 b 之間的距離，可以用 (Xa-Xb)平方加(Ya-Yb)平方再開個根號來計算。好像世界就該是這樣運轉的，對嗎?

其實不是每個空間都長這樣的。們在前面說過兩點之間的距離有很多種計算方式，但是不是每個距離都是這樣的。去看看微生物生態裡最常用到的 Bray-Curtis dissimilarity 或 UniFrac distance 你就會發現，它們之間的關係不符合 Euclidian distance。拿 Bray-Curtis dissimilarity 來看好了。假設 A 樣本在 OTU1, OTU2, OTU3, OTU4, OTU5 分別有 100, 100, 100, 100, 100 隻(總共 500 隻)，B 樣本有 0, 250, 150, 100, 50 隻(總共 550 隻)。在 OTU1, OTU2, OTU3, OTU4, OTU5 分別相差 100, 150, 50, 0, 100 隻。計算出來的分子是 100+150+50+0+100，分母是 500+550，計算出來的距離是 0.33。你可以自己玩玩看。而 A 樣本跟原點的距離是 1 (自己算)，B 樣本跟原點的距離也是 1，A-B-原點是個 1, 1, 0.3 的三角形。如果算歐幾里德距離，三邊距離分別是 224, 312, 212 的三角形，形狀完全不一樣了嘛。

這有什麼問題呢? 我們多變值統計方法裡，ordination 是個常用的呈現結果的方法。這個方法把所有樣本放在同一個圖裡，每個樣本的組成用一個點來代表，然後讓你用在圖上各點彼此是否靠近來判斷哪個樣本的菌相組成跟另一個樣本的菌相組成很接近。是個幫助你理解這個有很多變因的數據的好方法，因為看兩個點是否接近，跟比較兩個樣本裡每一個 OTU 的數量誰多誰少，要簡單而容易理解得多。

不幸的是，你容易理解的 X-Y 平面圖，是個歐幾里德平面，上面每個點之間的距離都是以歐幾里德距離來呈現的。而你原本估算距離用的指數，又因為生物學上的理由而不能用歐幾里德距離。有沒有大悲劇了?

還好數學上有辦法解決。常用方法裡有兩招，一個是 principal coordinate analysis (PCoA)，一個是 non-metric multidimensional scaling (nMDS)。這兩種方法都是需要你輸入各點間的距離當原料來分析，所以不管是怎麼算出來的距離都可以餵給它分析。兩種方法不同的地方是，PCoA 是擷取原來數據裡的符合歐幾里德距離的部份下來分析，只看一部份資源當然不足，但是至少有一部份能用，或許可以看出整體趨勢。而 nMDS 方法則是不管原始數據裡的數值，全部先依數值大小排名，再用排名當數據來分析(第 50 名就是 50)，這樣就符合歐幾里德距離，解決問題了。這樣做的好處是所有數據都用得到，壞處是如果 A1=100 A2=5 A3=4，你只能保留 A1>A2>A3 的次序，A1 遠大於另外兩個的數值資訊就沒有了。