RでMDS（多次元尺度構成法）②非計量MDS（NMDS）

veganパッケージのmetaMDS()関数を使ってNMDSを実行します。NMDSは、ランダムな初期配置から計算を始めるため、set.seed()で乱数を固定しておきます。

# 乱数を固定して結果を再現可能にする
set.seed(123)

# varespec（種構成データ）のみを使ってNMDSを実行
# 距離はブレイ・カーティス非類似度（デフォルト）
nmds_result <- metaMDS(varespec)

# 実行中の出力例
Square root transformation
Wisconsin double standardization
Run 0 stress 0.1843196 
Run 1 stress 0.2245479 
Run 2 stress 0.215148 
Run 3 stress 0.2324692 
Run 4 stress 0.1985581 
Run 5 stress 0.1843196 
... New best solution
... Procrustes: rmse 3.634903e-05  max resid 0.0001465453 
... Similar to previous best
Run 6 stress 0.1974408 
Run 7 stress 0.3840757 
Run 8 stress 0.2234787 
Run 9 stress 0.1843196 
... New best solution
... Procrustes: rmse 1.889957e-05  max resid 7.643299e-05 
... Similar to previous best
Run 10 stress 0.1852397 
Run 11 stress 0.2090733 
Run 12 stress 0.2066172 
Run 13 stress 0.1982376 
Run 14 stress 0.196245 
Run 15 stress 0.18584 
Run 16 stress 0.2085949 
Run 17 stress 0.1985582 
Run 18 stress 0.1985581 
Run 19 stress 0.2079056 
Run 20 stress 0.2126569 
*** Best solution repeated 1 times

metaMDS()は、最適な解を見つけるために異なる初期値から計算を繰り返します(今回はRun0~20の21回。デフォルトで20回以上)。

# 結果の概要を表示
nmds_result 

# 出力結果
Call:
metaMDS(comm = varespec) 

global Multidimensional Scaling using monoMDS

Data:     wisconsin(sqrt(varespec)) 
Distance: bray 

Dimensions: 2 　#1
Stress:     0.1843196 　#2
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries　#3
The best solution was from try 9 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘wisconsin(sqrt(varespec))’ 

#1 結果は2次元マップとして作成されました。
#2 ストレス値は0.184で、0.2を下回るため、解釈に使える良いモデルです（※）。
#3 21回の試行の中で、最も低いStress値を持つ安定した解が1種類見つかりました。

※💡ストレス値とは？
NMDSでは、元のデータの距離の大小関係を2次元あるいは3次元のマップでどれだけ忠実に再現できたか、その当てはまりの悪さを示す指標としてストレス値が計算されます。
この値は当てはまりの良さを示す数値として用いられ、0に近いほど良いモデルです。以下のような目安で判断されます※。
0.025: 素晴らしい
0.05: よい
0.1: まずまず
0.2: 悪い

※参考図書　金 明哲, Rによるデータサイエンス
※参考文献　土居 秀幸, 岡村 寛, 生物群集解析のための類似度とその応用 : Rを使った類似度の算出、グラフ化、検定, 日本生態学会誌, 2011, 61 巻, 1 号, p. 3-20,  https://doi.org/10.18960/seitai.61.1_3,

今回のストレス値は0.2を下回っており、信頼できる結果だと判断できます。

4-1. ステップ①：サンプル（地点）のみをプロット

サンプル（地点）の関係性だけを見て、全体像を掴みます。

# 地点のみをプロット（テキストラベルで表示）
plot(nmds_result, display = "sites", type = "text") # nmdsプロットをサンプル地点の番号としてを記載

このマップ上で距離が近い地点どうしは、下層植物の構成が似ていることを意味します。例えば、中央付近の地点19，23は互いに似ていると考えられます。

4-2. ステップ②：サンプルと種を追加してプロット

次に、地点の配置にどの種が影響しているのかを見るため、種の情報も追加します。

# 地点（黒）と種（赤）を両方プロット
plot(nmds_result, display = c("sites", "species"), type = "text")

地点の近くにある種は、その地点の植生を特徴づける種であると解釈できます。例えば、左上の地点2の近くにはCladstelがあり、これらの種がこのグループの植生に強く関連していることが考えられます。

# envfit()で環境変数とNMDS配置の相関を計算
set.seed(123) # 結果の再現性を確保
env_fit_result <- envfit(nmds_result ~ ., data = varechem)
env_fit_result

# 出力結果
***VECTORS

            NMDS1    NMDS2     r2 Pr(>r)    
N        -0.05039 -0.99873 0.2081  0.073 .  
P         0.68702  0.72664 0.1755  0.135    
K         0.82729  0.56178 0.1657  0.158    
Ca        0.75015  0.66127 0.2811  0.037 *  
Mg        0.69675  0.71731 0.3494  0.013 *  
S         0.27626  0.96108 0.1774  0.133    
Al       -0.83776  0.54603 0.5155  0.001 ***
Fe       -0.86197  0.50695 0.4001  0.007 ** 
Mn        0.80233 -0.59688 0.5322  0.002 ** 
Zn        0.66518  0.74668 0.1779  0.128    
Mo       -0.84876  0.52878 0.0517  0.569    
Baresoil  0.87210 -0.48933 0.2494  0.051 .  
Humdepth  0.92636 -0.37664 0.5590  0.001 ***
pH       -0.79908  0.60123 0.2624  0.044 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation: free
Number of permutations: 999

# プロットを描画
# 地点（黒）と種（赤）のマップを作成
plot(nmds_result, display = c("sites", "species"), type = "text")

# 環境変数の矢印（青）を重ね描き
plot(env_fit_result, col = "blue", p.max=0.05)

# 地点の座標を抽出
site_scores <- as.data.frame(scores(nmds_result, display = "sites"))
site_scores$site_name <- rownames(site_scores)

# 種の座標を抽出
species_scores <- as.data.frame(scores(nmds_result, display = "species"))
species_scores$species_name <- rownames(species_scores)

# 環境変数の座標を抽出
vector_scores <- as.data.frame(scores(env_fit_result, display = "vectors"))
vector_scores$variable <- rownames(vector_scores)

# vector_scoresにp値の列を追加
vector_scores$p_value <- env_fit_result$vectors$pvals 
# p値が0.05未満の行だけを抽出して、新しいデータフレームを作成
sig_vector_scores <- vector_scores[vector_scores$p_value < 0.05, ]

library(ggplot2)
ggplot() +
 # グレーで種をプロット
 geom_text(data = species_scores, aes(x = NMDS1, y = NMDS2, label = species_name),
color = "gray60", size = 3) +
 # 地点（サンプル）を黒文字でプロット
 geom_text(data = site_scores, aes(x = NMDS1, y = NMDS2, label = site_name),
color = "black", size = 4) +
 # p<0.05の環境変数の矢印
 geom_segment(data = sig_vector_scores,
  aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
  arrow = arrow(length = unit(0.2, "cm")), # 矢印の先の矢じりをかき、大きさを指定
  color = "blue", linewidth = 0.7) + # 矢印の色と線の幅指定
 # p<0.05の環境変数のラベルをプロット（矢印の少し先に表示）
 geom_text(data = sig_vector_scores,
  aes(x = NMDS1 * 1.1, y = NMDS2 * 1.1, label = variable),
  color = "blue", size = 4) +
 theme_test() +
 coord_fixed() # 縦横比を1:1に固定して距離の歪みをなくす