1. トップページ
  2. Rによる群集組成の解析
  3. 群集の類別と指標種

Rによる群集組成の解析

群集の類別と指標種

クラスター分析を行い、得られたグループに特徴的な指標種を抽出します。

群集間の類似度を計算し、階層的クラスターを作成

あらかじめrequire(vegan)にてveganパッケージを読み込みます。veganには群集の多変量解析に便利な関数が含まれています。またrequire(cluster)にてclusterパッケージを読み込みます。clusterには類別のための関数が含まれています。

群集データをファイルから読み込みます。以下の例ではcsvファイルから読み込んでいます。header=TRUEを明示して種名を読み込んでいます。また、地点名をrow.names=1で指定します。

dat <- read.csv("spcdat.csv", header = TRUE, row.names=1)

次に、群集組成の違いを数値化するために、類似度指数を基にした距離行列を求めます。veganの関数vegdist()を用い、method="horn"を指定してMorisita-Horn非類似度、method="bray"を指定してBray-Curtis非類似度の距離行列をそれぞれ作成します。

dis.mh <- vegdist(dat, method="horn")
dis.bc <- vegdist(dat, method="bray")

得られた距離行列をもとに、clusterの関数hclust()を用いて階層的クラスターを作成します。クラスターの作成方法にはいくつかの方法がありますが、次の例ではmethod="ward.D2"を指定してWard法によるクラスターを作成します。Ward法は他の方法よりもまとまりのある複数のグループを作りやすいといわれています。

resmh.ward <- hclust(dis.mh, method="ward.D2")
plot(resmh.ward)

クラスター作成の方法に群平均法をつかたい場合は、method="average"を指定します。上の例では地点Ako1が外れ値なのですが、Ward法で無理やりグループに組み込まれてしまった感じがします。

この図(デンドログラム)を見て、全体を3グループに分割したいと思った場合、cutree()関数を使用し、グループ数をk=3で指定します。

resmh.ward.g <- cutree(resmh.ward, k=3)
resmh.ward.g
##    Ako1    Ako2    Ako3  Nisio1  Nisio2  Nisio3 Nakatu1 Nakatu2 Nakatu3 
##       1       2       2       1       1       2       3       1       1 
##  Banzu1  Banzu2  Banzu3    Uto1    Uto2    Uto3 
##       2       2       2       3       3       1

この場合、グループの分割数は主観的に決めました。より客観的にグループ数を決める方法としてシルエット値を求める方法があります。

asil <- numeric(nrow(dat))
for (i in 2:(nrow(dat)-1)){
  sil <- silhouette(cutree(resmh.ward, k=i), dis.mh)
  asil[i] <- mean(sil[,3])
}
asil
##  [1] 0.00000000 0.27436092 0.37874725 0.41602592 0.43325536 0.41073129
##  [7] 0.34779031 0.32437756 0.28972742 0.23109473 0.24977641 0.22718512
## [13] 0.09599773 0.06826518 0.00000000

ここで、平均のシルエット値が最大となるのは、グループ数5の時の0.4332554です。従って、5分割が良いとなりましたので、その分割とシルエット図を表示します。

resmh.ward.g <- cutree(resmh.ward, k=which.max(asil))
resmh.ward.g
##    Ako1    Ako2    Ako3  Nisio1  Nisio2  Nisio3 Nakatu1 Nakatu2 Nakatu3 
##       1       2       2       3       3       2       4       5       5 
##  Banzu1  Banzu2  Banzu3    Uto1    Uto2    Uto3 
##       2       2       2       4       4       3
sil <- silhouette(resmh.ward.g, dis.mh)
rownames(sil) <- row.names(dat)
plot(sil)

非階層的クラスターを作成する

群集の組成が階層的に類別できるとは限らないので、非階層的に適当な類別を探索するほうが良い場合もあります。

非階層的クラスターを作成する方法として、K-means法がよく知られています。これは変数間のユークリッド距離を仮定していますので、群集組成の非類似度から求めた距離とは合わない面があります。そこで、よりロバストだと考えられているK-medoids法を使います。これにはclusterの関数pam()を利用します。

asilmk <- numeric(nrow(dat))
for (i in 2:(nrow(dat)-1)){
  asilmk[i] <- pam(dis.mh, k=i, diss=TRUE)$silinfo$avg.width
}
asilmk
##  [1] 0.00000000 0.21662533 0.32562256 0.36650728 0.42029671 0.38719369
##  [7] 0.38006387 0.34541372 0.28972742 0.23109473 0.21786278 0.08558420
## [13] 0.09599773 0.06826518 0.00000000

ここで、平均のシルエット値が最大となるのは、グループ数5の時の0.4202967です。従って、5分割が良いとなりましたので、その分割とシルエット図を表示します。

resmh.pam <- pam(dis.mh, k=which.max(asilmk), diss=TRUE)
summary(resmh.pam)
## Medoids:
##      ID            
## [1,] "1"  "Ako1"   
## [2,] "10" "Banzu1" 
## [3,] "15" "Uto3"   
## [4,] "9"  "Nakatu3"
## [5,] "13" "Uto1"   
## Clustering vector:
##    Ako1    Ako2    Ako3  Nisio1  Nisio2  Nisio3 Nakatu1 Nakatu2 Nakatu3 
##       1       2       2       3       4       2       5       4       4 
##  Banzu1  Banzu2  Banzu3    Uto1    Uto2    Uto3 
##       2       2       2       5       5       3 
## Objective function:
##     build      swap 
## 0.2251717 0.2217902 
## 
## Numerical information per cluster:
##      size  max_diss    av_diss  diameter separation
## [1,]    1 0.0000000 0.00000000 0.0000000  0.7825552
## [2,]    6 0.4932491 0.21581136 0.7984151  0.5436638
## [3,]    2 0.6921108 0.34605540 0.6921108  0.6004724
## [4,]    3 0.5783309 0.36657508 0.9096986  0.5536461
## [5,]    3 0.1268453 0.08004966 0.1640254  0.5436638
## 
## Isolated clusters:
##  L-clusters: character(0)
##  L*-clusters: [1] 5
## 
## Silhouette plot information:
##         cluster neighbor   sil_width
## Ako1          1        5  0.00000000
## Banzu1        2        5  0.67233690
## Banzu3        2        5  0.62128995
## Ako2          2        5  0.56283193
## Banzu2        2        5  0.52444944
## Nisio3        2        4  0.39476346
## Ako3          2        3  0.35069105
## Nisio1        3        4  0.14466868
## Uto3          3        5  0.07093218
## Nakatu3       4        3  0.33196280
## Nakatu2       4        5  0.20581221
## Nisio2        4        3 -0.04255231
## Uto1          5        3  0.86179224
## Uto2          5        1  0.81415322
## Nakatu1       5        2  0.79131885
## Average silhouette width per cluster:
## [1] 0.0000000 0.5210605 0.1078004 0.1650742 0.8224214
## Average silhouette width of total data set:
## [1] 0.4202967
## 
## Available components:
## [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
## [6] "clusinfo"   "silinfo"    "diss"       "call"
plot(resmh.pam)