α 多様性・β 多様性

概要
1. データ
2. alpha多様性
- shannon
- simpson
- chao1
3. beta多様性
- Bray-curtis
- Morishita
- Horn
- jaccard
4. デンドログラム描画

概要

α多様性とβ多様性は、主に群集の多様性を測定・評価する際に使用される。α多様性は、ある特定の生息地やサンプル内での種の数や個体の均等性（均質性）を評価する。一方でβ多様性は、異なる生息地や群集間でどの程度種が共有されているか、または異なっているかを評価し、その類似性や差異を理解する場合に用いられる。
ここでは、それぞれの統計指標について、Rのbase関数を用いて計算する場合と、veganパッケージを用いて計算する方法について記述する。

参考) 日本生態学会誌 61：3 - 20（2011）

1. データ

(dat <- data.frame(A=c(100,90,80,20,10,0),B=c(150,100,50,0,0,0),C=c(0,0,0,50,100,150),D=c(10,9,8,2,1,0)))

##     A   B   C  D
## 1 100 150   0 10
## 2  90 100   0  9
## 3  80  50   0  8
## 4  20   0  50  2
## 5  10   0 100  1
## 6   0   0 150  0

2. alpha多様性

shannon

\(H' = -\sum p_i ln p_i\)

shannon <- function(x) { pi <- x/sum(x); -sum(pi * log(pi), na.rm = T) }
sapply(dat, shannon)

##        A        B        C        D 
## 1.373774 1.011404 1.011404 1.373774

# veganを使う場合
vegan::diversity(t(dat), index = "shannon")

##        A        B        C        D 
## 1.373774 1.011404 1.011404 1.373774

simpson

\(D=1-\sum p_i^2\)

simpson <- function(x){pi <- x/sum(x); 1-sum(pi^2) }
sapply(dat, simpson)

##         A         B         C         D 
## 0.7222222 0.6111111 0.6111111 0.7222222

# veganを使う場合
vegan::diversity(t(dat), index = "simpson")

##         A         B         C         D 
## 0.7222222 0.6111111 0.6111111 0.7222222

chao1

\(S_{obs} + F_1^2/2F_2\)
NGSデータの場合、シングルトン・ダブルトンは削除するのが一般的(エラー由来リードの可能性)

# データ全体でみてF1,F2がきまるのか、特定のコレクションないのF1,F2を計算するのか
## 全コレクションにおけるF1,F2
chao1 <- function(dat){
  chao <- function(x, idx1, idx2){
    obs <- sum(x!=0 & !is.na(x))
    N <- sum(x, na.rm = T)
    F1 <- ifelse(length(x[idx1])>1, x[idx1], 0)
    F2 <- ifelse(length(x[idx2])>1, x[idx2], 0)
    a <- (F1^2)/(2*F2) * (N-1)/N
    ifelse(is.na(a) | is.infinite(a), obs, obs + a )
  }
  idx1 <- apply(dat, 1, max) == 1
  idx2 <- apply(dat, 1, max) == 2
  sapply(dat, function(x){chao(x,idx1,idx2)})
}
## 各コレクションにおけるF1,F2
chao1 <- function(x){
  obs <- sum(x!=0 & !is.na(x))
  N <- sum(x, na.rm = T)
  S1 <- sum(x==1) ; S2 <- sum(x==2)
  a <- S1*(S1-1)/(2*(S2+1)) * (N-1)/N
  ifelse(is.na(a) | is.infinite(a), obs, obs + a )
}

sapply(dat, chao1)

## A B C D 
## 5 3 3 5

# veganを使う場合
vegan::estimateR(t(dat))

##                  A   B   C         D
## S.obs    5.0000000   3   3 5.0000000
## S.chao1  5.0000000   3   3 5.0000000
## se.chao1 0.0000000   0   0 0.2236068
## S.ACE    5.0000000 NaN NaN 5.4910000
## se.ACE   0.8944272 NaN NaN 0.6970226

3. beta多様性

Bray-curtis

bray <- function(x, y){ sum(abs(x-y))/(sum(x)+sum(y)) }
as.dist(outer(dat, dat, Vectorize(bray)))

##           A         B         C
## B 0.2000000                    
## C 0.9000000 1.0000000          
## D 0.8181818 0.8363636 0.9818182

# veganを使う場合
vegan::vegdist(t(dat), method = "bray")

##           A         B         C
## B 0.2000000                    
## C 0.9000000 1.0000000          
## D 0.8181818 0.8363636 0.9818182

Morishita

morisita <- function(x, y){
  lmdx <- sum(x*(x-1))/(sum(x)*(sum(x)-1))
  lmdy <- sum(y*(y-1))/(sum(y)*(sum(y)-1))
  m <- 1-2*sum(x*y)/((lmdx+lmdy)*sum(x)*sum(y))
  ifelse(m<0, 0, m)
}
as.dist(outer(dat, dat, Vectorize(morisita)))

##            A          B          C
## B 0.06038159                      
## C 0.93288440 1.00000000           
## D 0.00000000 0.02735012 0.93052501

# veganを使う場合
vegan::vegdist(t(dat), "morisita")

##            A          B          C
## B 0.06038159                      
## C 0.93288440 1.00000000           
## D 0.00000000 0.02735012 0.93052501

Horn

horn <- function(x, y){ 
  lmdx <- sum(x^2)/(sum(x)^2)
  lmdy <- sum(y^2)/(sum(y)^2)
  1-2*sum(x*y)/((lmdx+lmdy)*sum(x)*sum(y))
}
as.dist(outer(dat, dat, Vectorize(horn)))

##            A          B          C
## B 0.06666667                      
## C 0.93333333 1.00000000           
## D 0.00000000 0.06666667 0.93333333

# veganを使う場合
vegan::vegdist(t(dat),"horn")

##            A          B          C
## B 0.06666667                      
## C 0.93333333 1.00000000           
## D 0.00000000 0.06666667 0.93333333

jaccard

距離に変換

# 在/不在データの場合(バイナリデータに変換)
(bin_dat <- data.frame(lapply(dat, function(x){ifelse(x!=0, 1, 0)})))

##   A B C D
## 1 1 1 0 1
## 2 1 1 0 1
## 3 1 1 0 1
## 4 1 0 1 1
## 5 1 0 1 1
## 6 0 0 1 0

dist_jaccard <- function(x, y){
 1 - (sum(x == 1 & y == 1)/(sum(x) + sum(y) - sum(x == 1 & y == 1))) 
}
as.dist(outer(bin_dat, bin_dat, Vectorize(dist_jaccard)))

##           A         B         C
## B 0.4000000                    
## C 0.6666667 1.0000000          
## D 0.0000000 0.4000000 0.6666667

# veganを使う場合 オプションbinary = T　をつければデータ変換する必要なし
vegan::vegdist(t(dat), method = "jaccard", binary = T)

##           A         B         C
## B 0.4000000                    
## C 0.6666667 1.0000000          
## D 0.0000000 0.4000000 0.6666667

4. デンドログラム描画

par(mfcol=c(1,2))
# Bray-curtis
bray <- function(x, y){ sum(abs(x-y))/(sum(x)+sum(y)) }
dist_bray <- as.dist(outer(dat, dat, Vectorize(bray)))
# dist_bray <- vegan::vegdist(t(dat), method="horn")
plot(as.dendrogram(hclust(d = dist_bray, method = "ward.D2")))

# Horn
horn <- function(x, y){ 
  lmdx <- sum(x^2)/(sum(x)^2)
  lmdy <- sum(y^2)/(sum(y)^2)
  1-2*sum(x*y)/((lmdx+lmdy)*sum(x)*sum(y))
}
dist_horn <- as.dist(outer(dat, dat, Vectorize(horn)))
# dist_horn <- vegan::vegdist(t(dat), method="horn")

plot(as.dendrogram(hclust(d=dist_horn, method = "ward.D2")))