Rを使って植生調査のデータを分類する方法

植生調査のデータを組成表としてまとめ、Rを使って組みかえて、
植物群落と標徴種を見つける方法を紹介します。

※実際には、植物の性質を基に、手作業で再調整する必要があります。あくまで、その手助けということで。

データの準備

地点名と出現種と、その優占度（被度など）の一覧表を用意しておきます。（EXCELとかで作業）

  • 複数の層で出現する種は、その最大値、合計、平均とかを採用する。（自分で決めて下さい）
  • 被度を、r,+,1〜5で記録している場合は、rと+→0.1,1→2.5,2→15,3→37.5,4→62.5,5→87.5に変換する。（被度の百分率表示）

例　"組成表.csv"

st1 st2 st3 st4 st5 st6 st7 st8 st9 st10
ヒメヒラテンツキ 62.5 87.5
イネ 62.5 87.5
ダンドボロギク 37.5 15
ヒメヘビイチゴ 2.5 2.5 0.1
セイタカアワダチソウ 15 87.5 15 87.5 87.5

組成表の組みかえの手順

  1. データを読み込んで、MDSもしくはDCAで序列化します。（metaMDS,decorana）
  2. 序列化したスコアを基に、階層的クラスター分析で、種、地点ごとにクラスタリングします（hclust）
  3. クラスタリングされた結果を、2分割→3分割→、、、→希望の分割数　まで順番に行い、その度に、同じグループ順となるように地点を並び替えます。（cutree,order）
  4. 種も同様に並び替えます。
  5. クラスタ数ごとに標徴種（指標種）を計算します。（indval）
  6. 組みかえた結果と、標徴種の結果をCSVに書き出します。
    • 組み替えた結果は、同じグループに属する地点のgroup欄に同じ値が入っています。2分割→指定した分割数　のすべての場合の　グループが示されます。
    • 徴標種の結果は、indvalの値が大きい(1が最大)ほど、そのグループを特徴づけており、pvalの値が小さいほど、その結果の信頼性が高いことを示しています。

Rプログラム

はじめにveganとlabdsvのパッケージをインストールしてください。

### 分割数設定　任意の数に変更してください　###
gnum_st<-15 ##地点
gnum_sp<-10　##種

### データ読み込み　###
library(vegan) 
library(labdsv) #indval用
d<-read.csv("組成表.csv",row.names=1)
d[is.na(d)]<-0
d<-t(d)
col<-ncol(d)
row<-nrow(d)

### 序列化 ###

#MDSの場合
d.mds<-metaMDS(d,zerodist="add") ##### zerodist="add" を追加

# 有る無しMDSの場合
#d[which(d>0,arr.ind=TRUE)]<-1 　
#d.mds<-metaMDS(d,dist="jaccard")

# DCAの場合
#d.mds<-decorana(d)

### クラスタリング ###
sco<-scores(d.mds,display="sites") #地点の序列化スコア
clus<-hclust(dist(sco),"average") #群平均法で

plot(clus) #グラフ表示
rect.hclust(clus,gnum_st,)

### グループ分けした序列を　グラフ表示 ###
windows()
k<-cutree(clus,gnum_st)
          #####　以下の1行プログラムミスがあったので修正 #####
ordipointlabel(d.mds,display="sites",col = c(2:(gnum_st+2))[k],pch = c(1:gnum_st+1)[k])
ordiellipse(d.mds,k, display="sites",kind="sd", conf=0.9,lwd=1,lty=2, col="black")

### 地点の分割数ごとにリストの並び替えと指標種を計算 ###
#変数初期化
li<-list(NULL)
me<-list(NULL)
ordtxt<-NULL
d2<-d #出力のためのデータ

#最大数に分割したときの順番を取得
tmp<-cbind(clus$order,c(1:row))
tmp<-tmp[order(tmp[,1]),]
tmp<-tmp[,2]

#分割数分だけループ
for( i in 1:gnum_st) {
 ifelse(i==gnum_st,me[[i]]<-tmp,me[[i]]<-cutree(clus,i+1))
 ordtxt[i]<-paste("me[[",i,"]]")#orderのための文字列
　d2<-cbind(d2,me[[i]])
 colnames(d2)[col+i]<-paste("group",i+1,sep="")

 dul<-indval(d,me[[i]]) #指標種分析 
 dul<-data.frame(numeric.sp=names(dul$maxcls), community=dul$maxcls, indval=round(dul$indcls,3), pvalue=round(dul$pval,4))
 dul<- dul[order(dul$community, dul$indval, dul$pvalue, decreasing=T),][,-1]
 li[[i]]<-dul #指標種分析結果
}
ord<-eval(parse(text=paste("order(",paste(ordtxt,collapse=","),")")))
d2<-d2[ord,]#並び替え

d2<-t(d2) ###地点、種入れ替え

### 種も同様に並び替え　###
sco<-scores(d.mds,display="species")
clus<-hclust(dist(sco),"average")
me<-list(NULL)
ordtxt<-NULL
tmp<-cbind(clus$order,c(1:col)) #クラスターの順番を表の並び替えに利用するため
tmp<-tmp[order(tmp[,1]),]
tmp<-tmp[,2]

for( i in 1:gnum_sp) {
 ifelse(i==gnum_sp,me[[i]]<-tmp,me[[i]]<-cutree(clus,i+1))
 ordtxt[i]<-paste("me[[",i,"]]")#orderのための文字列
　d2<-cbind(d2,me[[i]])
 colnames(d2)[row+i]<-paste("group",i+1,sep="")
}

# 並び替えの前に　必要ない行と列にNA代入
for(i in 1:gnum_st){
 for(j in 1:gnum_sp){
  d2[col+i,row+j]<-NA
 }
}
ord<-eval(parse(text=paste("order(",paste(ordtxt,collapse=","),")")))
d2<-rbind(d2[ord,],d2[(col+1):(col+gnum_st),])#並び替え

d2<-d2[,-(row+gnum_sp)]　#並び替えのための最終列を削除
d2<-d2[-(col+gnum_st),] #並び替えのための最終行を削除

### 並び替えた結果と指標種を書き出し
write.table(d2,"結果.csv",sep=",",row.names=T,col.names=T)
for(i in 1:gnum_st) write.table(li[[i]],paste("indval",i,".csv",sep=""),sep=",",row.names=T,col.names=T)

メモ（自分用）

  • metaMDSとdecoranaで、データによっては結果がまあまあ違う。比較してみた方がいいかも。
  • metaMDSのnoshareのオプションでtoo long のエラーを制御できるけど、意味はまだ理解していない。too longのwarningがでるデータだと結果があまりよくない？
  • 植生データベースを使って反復平均法と結果を比較しようと試したけど、CSVを上手く登録できなかったので挫折。
  • winTWINSを試そうと思ったけど、データの入力形式が面倒なので、挫折。
  • 指標種分析の代わりに、決定木を使ったらどうだろう？基準変数はグループ番号で、説明変数を優占度にする。
  • 序列化に自己組織化マップを使ったらどうだろう？

参考文献