生物系学生のメモ帳

Rで解析 クラスター解析と指標種分析まで

ここでは、ワォード法によるクラスター解析と、Indval法による指標種分析を行うコード例を紹介します。解析したいデータに合わせて、適宜変更しながら使ってみてください（エラーの理由とか対処法とかはChatGPTなんかを使いながら対応していくのがオススメ）。

データセットは、以下の様なものを想定しています。（データセット例は超適当なので気にしないで下さい）

（データセット例）

標準和名＼地点名　　沖縄　　和歌山　　東京　　宮城　

マルタ　　　　　　　　　　　　　　　　　　　　5

シロザケ　　　　　　　　　　　　　　　　　　　10

ボラ　　　　　　　　15　　　15　　　　50　　　30

コボラ　　　　　　　10　　　　　　　　　　　　　

メナダ　　　　　　　　　　　20　　　　15　　　20

クロダイ　　　　　　　　　　50　　　　30　　　30

ミナミクロダイ　　　50　　　　　　　　　　　　　

リボンスズメダイ　　10　　　　　　　　　　　　　

カワスズメ　　　　　70　　　　　　　　　　　　　　　

ゴマフエダイ　　　　30　　　10　　　　　　　　　　

ゴマアイゴ　　　　　30　　　　　　　　　　　　　　　

オキナワフグ　　　　20　　　　　　　　　　　　　　

あんまりスマートなコードではないかもしれません。「こうした方が良いよ！」っていうのを思いついた方は、コメント等で指摘していただけるととてもありがたいです。

以下、コード例。

library(tidyverse)

library(vegan)

library(cluster)

# 作業ディレクトリを指定（ファイルが保存されているフォルダ）

setwd("作業したいファイルを扱うフォルダのディレクトリ")

# ファイルの読み込み

Result = read.csv("扱いたいファイル.csv", fileEncoding = "CP932")

DATASET <- Result

# NAを0に変換

DATASET2 = DATASET

DATASET2[is.na(DATASET)] = 0

# クラスタリングに使用する数値データを作成（標準和名列を除外）

numeric_data <- DATASET2[, -1]

# 各列の値を対数変換（Log10(X + 1））をDATASETに適用

DATASET3 = numeric_data

columns_to_transform = 1:26

DATASET3[, columns_to_transform] = log10(DATASET3[, columns_to_transform] + 1)

DATASET4 = DATASET3

# データフレームの行と列を入れ替える

DATASET4 <- t(DATASET4)

DATASET5 = DATASET4

＃ここから、明示するまでは大半の人が必要ない操作だと思う

#20行目を削除（何も取れなかった時のデータのため、空データがあると解析出来ないみたい）

transposed_data <- DATASET5[-20, ]

#72列目を削除（何か勝手に追加されてた空の列。これも、空データがあると解析できないみたいなので削除しておく）

transposed_data2 <- transposed_data[,-72 ]

＃大半の人が必要ないと思われる操作はここまで。

# 以下、クラスタリング

dis.mh = vegdist(transposed_data2, method="horn")

dis.bc = vegdist(transposed_data2, method="bray")

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

plot(resmh.ward)

#クラスター数の最適解を推定

asil = numeric(nrow(transposed_data2))

for (i in 2:(nrow(transposed_data2)-1)){

sil <- silhouette(cutree(resmh.ward, k=i), dis.mh)

asil[i] <- mean(sil[,3])

}

asil

resmh.ward.g <- cutree(resmh.ward, k=which.max(asil))

resmh.ward.g

sil <- silhouette(resmh.ward.g, dis.mh)

rownames(sil) <- row.names(transposed_data2)

plot(sil)

#出力

png("silhouette_plot_Log.png", width = 1500, height = 1500, res = 200)

plot(sil)

dev.off()

# シルエット値が最も高かったクラスター数を取得

optimal_clusters <- which.max(asil)

print(paste("最適なクラスター数:", optimal_clusters))

# デンドログラムを再プロット

plot(resmh.ward, hang = -1, cex = 1, lwd = 2)

#出力

# PNG形式で保存

png("保存したいファイル名.png", width = 10, height = 13, units = "in", res = 1200)

# マージンを調整

par(mar = c(10, 4, 2, 2))

# 正しいオブジェクトでプロット

plot(resmh.ward, hang = -1, cex = 2, lwd = 3)

# ファイルを閉じる

dev.off()

# JPEG形式で保存

jpeg("保存したいファイル名.jpg", width = 600, height = 750)

# 正しいオブジェクトでデンドログラムをプロット

plot(resmh.ward, hang = -1, cex = 5, lwd = 3)

# ファイルを閉じる

dev.off()

#グループの指標種をIndVal法で抽出する

# labdsvパッケージをインストール（初回のみ）

install.packages("labdsv")

# labdsvパッケージを読み込む

library(labdsv)

# iva <- indval(dat, resmh.pam$clustering, numitr = 5000)の呼び出し

iva <- indval(transposed_data2, resmh.ward.g, 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(transposed_data>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),])

write.csv(fdig, "C保存したい場所のディレクトリ\\保存したいファイル名.csv",row.names = TRUE)

#ここまで。何か面白い発見がありますように。

指標種に関しては、ここでは行数で出てくる（解析の時に数値データのみのファイルにするために和名列を削除したため）。行数と和名列を照らし合わせてみて、指標種が何なのかを調べることになります。