1. トップページ
  2. Rによる群集組成の解析
  3. 群集組成による配置図:間接傾度分析

Rによる群集組成の解析

群集組成による配置図:間接傾度分析

群集組成の良く似た地点同士は近くに、似ていない地点は遠くになるように、空間上に地点を配置します。図示すると地点間の関係が認識しやすくなります。さらに、その配置図に対する環境変数の相関を求めます。

群集間の類似度を計算し、nMDSプロットを作成

nMDSは非計量多次元尺度法と呼ばれます。多次元で表された多数の点を、お互いに類似した点同士は近く、類似していない点同士は遠くなるように、なるべく低い次元(二次元平面など)に配置し直す方法です。群集データの場合は、地点間の類似性は類似度指数で表します。

あらかじめrequire(vegan)にてveganパッケージを読み込みます。veganには群集の多変量解析に便利な関数が含まれています。
群集データをファイルから読み込みます。以下の例ではcsvファイルから読み込んでいます。header=TRUEを明示して種名を読み込んでいます。また、地点名を行の名前にするためにrow.names=1を指定します。

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

次に、metaMDS()で類似度指数を指定して、nMDSの計算を行い、結果を図示します。
以下の例ではMorisita-Horn類似度を用いています。Morisita-Horn類似度指数は、木元Cπ指数と同等の指数です。

res1.mds <- metaMDS(dat, distance="horn", autotransform = FALSE, trace = 0)
plot(res1.mds, type="t", display="sites")

すると、nMDSの配置図が表示されます。

元々は種類数の分だけ次元があるので102次元のデータです。2次元の平面に配置するには相当無理をしていると考えられます。どの程度無理をしているかを表す指標がstress値です。計算の出力を見てみます。

res1.mds
## 
## Call:
## metaMDS(comm = dat, distance = "horn", autotransform = FALSE,      trace = 0) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     dat 
## Distance: horn 
## 
## Dimensions: 2 
## Stress:     0.1590454 
## Stress type 1, weak ties
## Two convergent solutions found after 20 tries
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'dat'

この場合、stress値は0.1590454となりました。一般的にstress値が0.2以上だと良くないといわれています。その場合は次元を上げてnMDSをおこないます。三次元にするには引数としてk=3を指定し、metaMDS(dat, distance="horn", k=3, autotransform = FALSE, trace = 0)とすれば良いです。このデータで三次元のnMDSを行うと、stress値は0.0948034と改善されます。

群集間の類似度を計算し、PCoAプロットを作成

上のnMDSの例では、群集間の類似度の計算と図示を一つの関数metaMDSで行いました。とても便利な方法です。次の例では、サンプル間の群集組成の違いを距離(非類似度)で表した行列を作成し、その距離行列から二次元のプロットを作成します。類似度は、逆数をとる、もしくは1から引く等の操作をすると非類似度となり、指数の性質により適当な操作をおこなえば「距離」として使うことができます。

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

上の行でMorishita-Horn非類似度、下の行でBray-Curtis非類似度の距離行列をそれぞれ作成しました。ここで、metaMDSを使うと先の例と同様のnMDSプロットが作成できます。次の例ではBray-Curtis非類似度の距離行列を使ってnMDSの二次元プロットを作成します。metaMDSは生データと距離行列のどちらからでもnMDSプロットを作成することができるのです。

res2.mds <- metaMDS(dis.bc, trace=0)
plot(res2.mds, type="t", display="sites")

では、PCoAプロットを作成します。cmdscale関数を使います。この関数は、距離行列を引数に使用します。次の例では引数にdis.mhを使ってMorishita-Horn非類似度によるプロットを作成しました。Bray-Curtis非類似度を使いたければ引数に上で計算したdis.bcを使ってください。k=2で第2軸までの計算を指定しています。

res.pco <- cmdscale(dis.mh, k=2, eig=TRUE)
ordiplot(res.pco, type="t")

PCoAは主座標分析や計量多次元尺度法と呼ばれ、後で解説するdbRDA (distance-based redundancy anaysis)にて行う環境傾度分析の基礎になります。

では、nMDSで作成した配置図とPCoAで作成した配置図はどの程度異なっているのでしょうか。配置図の比較はprocrustesテストprotest()で行います。

pro <- protest(res1.mds, res.pco)
pro
## 
## Call:
## protest(X = res1.mds, Y = res.pco) 
## 
## Procrustes Sum of Squares (m12 squared):        0.2582 
## Correlation in a symmetric Procrustes rotation: 0.8613 
## Significance:  0.001 
## 
## Permutation: free
## Number of permutations: 999
plot(pro)

Procrustes testでは有意確率0.001で異なっているという結果になりました。

群集配置図と環境変数の関連を図示

群集の配置図では、近くに配置された地点は良く似た群集組成、遠くに配置された地点はあまり似ていない群集組成を表しています。この群集組成の違いが、環境変数と何らかの関連があるかどうか、調べてみます。

まず、環境変数のデータを読み込みます。

edat <- read.csv("envdat.csv", header=TRUE)

次に、nMDSの配置図との相関を調べます。関数envfit()を使います。

ef <- envfit(res2.mds, edat, permu=4999)
ef
## 
## ***VECTORS
## 
##         NMDS1    NMDS2     r2 Pr(>r)  
## PSU   0.29621 -0.95512 0.1732 0.3232  
## PH    0.86206 -0.50680 0.4043 0.0454 *
## ORP3  0.25097  0.96799 0.2325 0.2066  
## ORP5 -0.01330  0.99991 0.3348 0.0926 .
## ORP8 -0.00882  0.99996 0.3121 0.1120  
## IL   -0.22412  0.97456 0.3402 0.0886 .
## Chla -0.67474  0.73806 0.0363 0.7982  
## MdG  -0.72298 -0.69087 0.3608 0.0348 *
## PM    0.16020 -0.98708 0.2993 0.1084  
## TN   -0.04183  0.99912 0.1270 0.4814  
## TOC  -0.91193 -0.41035 0.1318 0.4850  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 4999
ordiplot(res2.mds, type="t", display="sites")
plot(ef, p.max=0.1)

先の解析で環境変数のうち、PH, ORP5, IL, MdGの相関が高いということがわかりました。配置図上でこれら変数の増加する方向をベクトル表示しました。次の例では、これら4個の変数の等高線表示を行います。

sgev <- names(ef$vectors$r[ef$vectors$pvals<0.1])
ordiplot(res2.mds, display="sites")
for (i in 1:length(sgev)){
 fmi <- formula(paste("res2.mds~", sgev[i]))
 tmp <- with(edat, ordisurf(fmi, add=TRUE, col=i+1))
}

等高線表示により環境変数と群集組成の配置図との関係が、ベクトル表示よりも詳細にわかります。複数の等高線が重なり合ってに見にくい時は、with(edat, ordisurf(res2.mds, Ph, add=TRUE, col="red"))のように、変数ごとに作図したほうが良いです。

等高線の間隔がふぞろいで、しかも曲線を描いていることから、有意な環境変数であっても、配置図に対して環境変数が単純な線型の関係ではないことがわかります。また、山形の等高線を描く環境変数では、線型ではない深い関係がある可能性を示しています。


参考文献

  • Borcard D, Gille F, Legendre P (2011) Numerical Ecology with R. Springer
  • Legendre P, Legendre L (2012) Numerical Ecology. Elsevier
  • Magurran AE, McGill BJ (2011) Biological Diversity. Oxford Univ.

Rによる群集組成の解析に戻る