HackMD比較物種豐富度的手段：稀釋曲線與樣本涵蓋率估計
tags: 生態學 學習筆記
在生物多樣性的估計量整理（上）：緒論與 α 多樣性一文中，提到期望物種數的概念，是為了比較不同的群集物種豐度而發展出來的。因為不同群集、或者不同調查趟次所得的總個體數不同，又或者努力量不同，就算估計出了有效物種數，也沒有辦法直接比較生物多樣性的高低，因此需要標準化之後，方可進行比較。
美國海洋無脊椎動物學家 Howard L. Sanders 在 1968 年一篇比較不同海底無脊椎動物群集豐度的文章中，提出了稀釋化（rarefaction）的作法[1]，透過比較各樣點採獲的物種個體數與調查到物種百分比之間的關係，呈現其物種多樣性的情況。
其實在 Sanders 之前，已經有許多學者，企圖以計量方法呈現調查努力量（通常以調查到的個體數來代表）與調查到的物種數之間的關係。下表取自 Bobrowsky and Ball 1989，稍微調整代數符號：
| 計量方法 | 來源 |
|---|---|
| \(\hat{S}_{1}=\alpha\left(\ln\left(1+\frac{n}{\alpha}\right)\right)\) | Fisher, Corbet, and Williams (1943) [2] |
| \(\hat{S}_{2}=y_{0}\hat{\sigma} \sqrt{2\pi}\) | Preston (1948) [3] |
| \(\hat{S}_{3}=2.07\left(\frac{n}{X_{r}}\right)^{0.262}\) | Preston (1962a) [4], Preston (1962b) [5] |
| \(\hat{S}_{4}=2.07\left(\frac{n}{X_{r}}\right)^{0.262}A^{0.262}\) | Macarthur (1965) [6] |
| \(\hat{S}_{5}=kA^{d}\) | Kilburn (1966) [7] |
| \(\hat{S}_{6}=\frac{an}{1+bn}\) | de Caprariis, Lindemann, and Collins (1976) [8] |
| \(\hat{S}_{7}=\sum_{i=1}^{S_{obs}}1-\frac{\binom{n-X_i}{m}}{\binom{n}{m}}\) | Hurlbert (1971) [9] |
\(S\)、\(S_{obs}\)、\(n\)、\(m\)、\(X_i\) 之涵義見下段；\(\alpha\) 代表 Fisher 斜率常數、\(y_0\) 代表模型中的物種或類型數量、\(X_{r}\) 代表群集中最稀有的物種個體數、\(A\) 代表調查的面積、\(\hat{\sigma}\) 代表估計的標準差、\(k, d, a, b\) 為迴歸分析中會得到的係數。
經過其他生態學家與統計學家的努力，對於稀釋化的方法已經漸漸產生了標準的統計模型。
稀釋化方法的模型
假設一個群集已知有 \(S\) 個物種，各物種的相對豐富度為 \(p_i\)，以歸還抽樣法從群集中隨機取 \(n\) 個體，令 \(X_i\) 是第 \(i\) 個物種在樣本中出現的次數，則可以假設 \(X_1,X_2,...,X_S\) 的分配為：
\[P(X_1,X_2,...,X_S)=\frac{n!}{x_1,x_2!,...,x_S!}p_1^{x_1}p_2^{x_2}...p_S^{x_S}\]
當抽樣數為 \(m\) 時，由 Good (1953)[10]推導出的定理，可以得此時看到的物種數期望值 \(E[S_{m}]=S-\displaystyle\sum_{i=1}^{S}(1-p_i)^m\)。
對於生態調查資料，假設調查到的個體數為 \(n\)、調查到的物種數為 \(S_{obs}\)，Smith and Grassle (1977)[11] 由 Hurlbert (1971) 的式子，證明在資料中隨機抽樣數為 \(m\ (m<n)\) 時，觀察到此時物種數 \(S_{m}\) 的最小變異數不偏估計量如同下式：
\[\hat{S}_{m}=S_{obs}-\displaystyle\sum_{X_i>0} \frac{\binom{n-X_i}{m}}{\binom{n}{m}}\]
只要把調查資料放入式中，就可以估計調查的個體數稀釋到 \(m\) 時，物種數的期望值為何；如果不同的調查個體數不相同，只要用此方法把期望物種數「稀釋」至相同，就可以比較不同群集，或者相同樣點、不同時間之間的調查物種豐度高低。把對任意 \(m<n\) 的物種數估計對 \(m\) 做圖，可以得到一條稀釋曲線，任何對 \(m<n\) 的物種數估計都是對稀釋曲線做內插法。
稀釋化方法的外推模型
有了稀釋化的方法，我們在比較多個群集的物種豐度時，就可以藉由稀釋曲線的內插，將各群集的個體數調整到與規模最小的群集相同來進行比較。從另一個角度來看，這樣就捨去了個體數比較多的那些群集多出來的物種；這些資料如果能再拿來利用，那該有多好。
在稀釋曲線取樣模型的基礎下，沈宗荏博士等人在 2003 年發展出了稀釋曲線的進一步取樣統計模型[12]，考慮在構成稀釋曲線的取樣下，估計已經取完全部的 \(n\) 個體後，若再取得 \(m^*\) 個體時，可以看到多少物種數。利用 Good (1953) 與 Smith and Grassle (1977) 的研究結果，可以得到此時物種數不偏估計量為：
\[\hat{S}_{n+m^*}=S-\displaystyle\sum_{i=1}^{S}(1-p_i)^{n+m^*}\]
不過如同抽樣數為 \(m\) 的模型一樣，這裡的估計量並沒有辦法估計變異數。
估計稀釋化模型的變異數
要估計稀釋曲線外插的變異數，事實上就是需要把估計取樣中沒有出現，亦即出現次數為 0 的物種 \(f_0\) 考慮進去。其實 Shen(2003) 已經有提出，後來經過趙蓮菊老師修飾之後[13]表示如下：
\[\hat{S}_{(n+m^*)}=S_{obs}+\hat{f_0}\left[1-\left(1-\frac{f_1}{n\hat{f_0}+f_1}\right)^{m^*}\right]\]
要估計取樣中沒有出現的物種數，聽起來是個不可能的任務，不過 Colwell 等人 (2012)[14] 認為可以使用們在生物多樣性的估計量整理（上）：緒論與 α 多樣性一文中提到的 Chao1 豐富度指數估計方法，其中就有使用出現 1 與 2 次的物種數來估計 \(f_0\) 的下界：
\[\begin{split}
\hat{f_0}&=\frac{n-1}{n}\frac{f_1^2}{2f_2}\text{,}\ f_2>0\\
\\
\text{or }\ \hat{f_0}&=\frac{n-1}{n}\frac{f_1(f_1-1)}{2(f_2+1)}\text{,}\ f_2=0
\end{split}\]
有了變異數估計量的下界，就可以配合拔靴法（bootstrap）定出外插範圍的 95% 信賴區間。趙蓮菊老師建議這邊的 \(m^*\leq n\)，意思是最多外插估計到調查個體數的兩倍就好，以免誤差過大。
透過標準化的稀釋方法估計樣本涵蓋率
一路看下來，我們知道稀釋曲線可以外推、可以透過拔靴法得到信賴區間，但若要用以比較基於不同個體數的調查樣本，是否有標準化的方法，以取得更為公允的估計？
在二戰期間，知名的密碼學家、計算機科學家圖靈（Turing）與其同事 Good 應破譯密碼的需求，已經發展了樣本涵蓋的概念。針對數量為 n 的一組樣本，其涵蓋度定義為：
\[C_{n}=\displaystyle\sum_{i=1}^{S}p_iI\ (X_i>0)\]
其中 \(p_i\) 為第 \(i\) 組成的比例，\(I\) 為樣本的指示函數（如充裕度或出現與否等等），\(X_i\) 為第 \(i\) 組成的出現次數。
以物種調查資料的角度，\(C_{(n)}\) 就是在樣本中看到的各物種機率總和。若要精確估計樣本完整性，Good 與 Turing 發現可以透過出現 1 次的物種比例來表示[15]：
\[\hat{C}_{n}=1-\frac{f_1}{n}\]
針對上述估計，趙蓮菊老師與 Jost (2012)[16] 引入出現 2 次的物種比例，提出更精確的估計為：
\[\hat{C}_{n}=1-\frac{f_1}{n}\left[\frac{(n-1)f_1}{(n-1)f_1+2f_2}\right]\]
如此一來，如同前文提到的物種數期望值，此時可推得對任何樣本數為 \(m\) 的樣本涵蓋 \(C_{(m)}\) 期望值 \(E[C_{(m)}]\) 為：
\[E[C_{m}]=1-\displaystyle\sum_{i=1}^{S}p_i(1-p_i)^m\]
接著，針對稀釋樣本數 \(m<n\) 的物種涵蓋 \(\hat{C}_{m}\) 最小變異數不偏估計量如下：
\[\hat{C}_{m}=1-\displaystyle\sum_{X_i\geq1} \frac{X_i}{n} \frac{\binom{n-X_i}{m}}{\binom{n-1}{m}},\ \ \ m<n\]
最後，若要多取樣 \(m^*\) 個物種進行外推，需要將未出現的物種一並考慮進去。考慮稀釋曲線的外推方法與前文提到的 \(\hat{f_0}\) 估計方法，最後可以得到樣本涵蓋率的稀釋曲線估計方法：
\[\hat{C}_{(n+m^*)}=1-\frac{f_1}{n}\left[\frac{(n-1)f_1}{(n-1)f_1+2f_2}\right]^{m^*+1}\]
當 n 很大的時候，此式會趨近於 \(\left(1-\hat{C_{n}}\right)exp\left[\frac{-2m^*f_2}{(nf_1)}\right]\)，最後可以得到稀釋曲線的期望物種數與期望樣本涵蓋率的關係為：\(1-\hat{C}_{n+m^*}=\hat{S}_{n+m^*+1}-\hat{S}_{n+m^*}\)。
整合期望物種數與樣本涵蓋率
若以期望物種數 \(\hat S\) 對樣本涵蓋率 \(\hat C\) 作圖，可以得到稀釋曲線與外插曲線，兩條曲線會交於 \((\hat C_n, S_{obs})\)，意即該抽樣／調查到的物種數（\(S_{obs}\)）的樣本涵蓋率（\(\hat C_n\)）。
因為樣本涵蓋率是經過標準化之後的數值，因此可以用以客觀比較不同樣本之間的物種豐富度。若將 \(m^*\) 取至建議的最大值（\(n\)），我們得以估計如果大幅增加調查努力量，擴大調查的個體數到 \(2n\) 時的物種數與涵蓋率期望值。
舉個例子：
以下圖 (c) 為例，在繪出樣區的樣本涵蓋率稀釋與外插曲線後，我們可以得知 old growth 與 second growth 兩樣區調查結果的樣本涵蓋率。透過稀釋曲線與外插法模型，可以估計出如果投入更多努力做調查，當調查到的個體數加倍時，old growth 樣區的樣本涵蓋率從 0.65 增加到 0.73，而 second growth 樣區的樣本涵蓋率則只增加了 0.03；結合期望樣本涵蓋率與期望物種數來看，可以比較兩個樣區在任意涵蓋率 <0.73 時，old growth 樣區的物種數內插值或期望值均比 second growth 樣區來得多，顯見 old growth 樣區具有較高的物種豐度。
| old growth | second growth | |
|---|---|---|
| 調查物種數 | 112 | 140 |
| 最大期望物種數* | 194 | 198 |
| 物種涵蓋率 | 0.65 | 0.93 |
| 最大期望涵蓋率* | 0.73 | 0.96 |
*即樣本數為取樣數 2 倍時的物種數與樣本涵蓋率期望值。
取自 Chao & Jost (2012)。
原本的稀釋插曲線，已經可以讓人很直觀地比較樣本間的物種豐富度；加上有限度的外插方法，可以把資料使用得更淋漓盡致。從上面的例子，我們可以思考，如果對不同特性的樣區以同樣努力量進行調查或規定調查到相同的個體數即停止，低多樣性的群集可以被充分調查，但高多樣性群聚的群集調查結果卻可能沒有代表性。
以樣本涵蓋率的估計方法，讓使用者得以將調查資料標準化，可以讓物種豐度的比較結果更有說服力，也可以據此彈性調整生態調查的努力量。
延伸閱讀
- 許皓捷(2019)。動物調查努力量應隨棲地與季節調整。台灣生物多樣性研究，21(2)，41-57。
- Gotelli, N. J., & Chao, A. (2013). Measuring and estimating species richness, species diversity, and biotic similarity from sampling data. In: Encyclopedia of Biodiversity, pp.195-211.
- 蘇黎世大學的課程網頁Measures of biodiversity。
🐕🦺2023.02.24
Sanders, H. L. (1968). Marine benthic diversity: a comparative study. The American Naturalist, 102(925), 243-282. ↩︎
Fisher, R. A., Corbet, A. S., & Williams, C. B. (1943). The relation between the number of species and the number of individuals in a random sample of an animal population. The Journal of Animal Ecology, 12(1), 42-58. ↩︎
Preston, F. W. (1948). The commonness, and rarity, of species. Ecology, 29(3), 254-283. ↩︎
Preston, F. W. (1962). The canonical distribution of commonness and rarity: Part I. Ecology, 43(2), 185-215. ↩︎
Preston, F. W. (1962). The canonical distribution of commonness and rarity: Part II. Ecology, 43(3), 410–432. ↩︎
MacArthur, R. H. (1965). Patterns of species diversity. Biological reviews, 40(4), 510-533. ↩︎
Kilburn, P. D. (1966). Analysis of the species‐area relation. Ecology, 47(5), 831-843. ↩︎
de Caprariis, P., Lindemann, R. H., & Collins, C. M. (1976). A method for determining optimum sample size in species diversity studies. Journal of the International Association for Mathematical Geology, 8(5), 575-581. ↩︎
Hurlbert, S. H. (1971). The nonconcept of species diversity: a critique and alternative parameters. Ecology, 52(4), 577-586. ↩︎
Good, I. J. (1953). The population frequencies of species and the estimation of population parameters. Biometrika, 40(3-4), 237-264. ↩︎
Smith, W., & Grassle, J. F. (1977). Sampling properties of a family of diversity measures. Biometrics, 33(2), 283-292. ↩︎
Shen, T. J., Chao, A., & Lin, C. F. (2003). Predicting the number of new species in further taxonomic sampling. Ecology, 84(3), 798-804. ↩︎
Chao, A., Colwell, R. K., Lin, C. W., & Gotelli, N. J. (2009). Sufficient sampling for asymptotic minimum species richness estimators. Ecology, 90(4), 1125-1133. ↩︎
Colwell, R. K., Chao, A., Gotelli, N. J., Lin, S. Y., Mao, C. X., Chazdon, R. L., & Longino, J. T. (2012). Models and estimators linking individual-based and sample-based rarefaction, extrapolation and comparison of assemblages. Journal of plant ecology, 5(1), 3-21. ↩︎
Good, I. J. (1953). The population frequencies of species and the estimation of population parameters. Biometrika, 40(3-4), 237-264. ↩︎
Chao, A., & Jost, L. (2012). Coverage‐based rarefaction and extrapolation: standardizing samples by completeness rather than size. Ecology, 93(12), 2533-2547. ↩︎
