- トップページ
- Rによる群集組成の解析
- 群集の類別と指標種
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 3sil <- 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)
グループの指標種をIndVal法で抽出する
群集が類別できると、それぞれのグループに特徴的に出現する種を特定できます。ここでは、IndVal法による指標種の抽出を行います。IndVal法では、グループごとに各種のIndVal値をもとめ、有意に高い値を示す種をそのグループの指標種とします。IndValの計算には、群集データと群集の類別結果が必要です。
ここでは、K-medoidsでの類別結果を用います。また、indval()関数はlabdsvパッケージにありますので、require(labdsv)で読み込んでおきます。
iva <- indval(dat, resmh.pam$clustering, numitr=5000)
gr <- iva$maxcls[iva$pval<=0.05]
iv <- iva$indcls[iva$pval<=0.05]
pv <- iva$pval[iva$pval<=0.05]
fr <- apply(dat>0,2,sum)[iva$pval<=0.05]
fdig <- data.frame(group=gr, indval=iv, pvalue=pv, freq=fr)
(fdig <- fdig[order(fdig$group, -fdig$indval),])## group indval pvalue freq
## Musculista_senhousia 1 0.9926342 0.0292 9
## Nectoneanthes_latipoda 1 0.9795918 0.0360 2
## Aonides_oxycephala 1 0.9354839 0.0370 6
## Spio_sp 3 0.8571429 0.0158 5
## Synidotea_sp 4 1.0000000 0.0040 3
## Mactra_veneriformis 5 0.8463950 0.0026 9各グループについて、有意(p<0.05)にindval値の大きい種が指標種として抽出されました。指標種のないグループもあります。
参考文献
- Borcard D, Gille F, Legendre P (2011) Numerical Ecology with R. Springer
- Legendre P, Legendre L (2012) Numerical Ecology. Elsevier
Rによる群集組成の解析に戻る