EcoPlate解析 (3)

５）多変量解析の基礎

ここではエコプレートのデータ（32次元のベクトル）をはじめ、群集の種組成など、生態学でよく出くわす多次元データ＝「多変量データ」に関する統計の基礎について解説します。ここでの解説の内容は、Tree Diversity Analysisというサイトの中の、Tree diversity analysis manualというマニュアルの第8章・10章の一部を参考にしています。より一貫した情報を得るにはこのマニュアルの熟読をお勧めします。

この講習では、先ほど作ったエコプレートのデータの一部(osaka_summary_integ_csv)を使って、多変量解析の「イメージ」を説明したいと思います。データをシンプルにしたりきれいな図を描いたりすることが可能な解析方法の背後には数学理論があり少なくとも線型代数の知識が必要なのですが行列計算や固有値など数学的説明はすべて省略します。

複雑なデータの可視化

では初めに、osaka_summary_integ.csvをたとえばエクセルで開いてみます（ついでにエクセル形式で保存しなおします）。このファイルは縦方向に行(column)が並び、横方向に列(row)が並んでいます。縦方向の各行V1-V31の31個の値が、各サンプルの31種類の炭素基質に対する代謝能力を示しています。問題は、サンプルごと31個の値全体の特徴を考慮して、サンプル間のばらつきをどうやって可視化するか、です。が、そんなの簡単です。エクセルで棒グラフ描けばよいですよね。上から下までの32個のサンプル全部を使ってグラフを描くと下の感じになります。

ぐちゃぐちゃしてますね。ちなみにmetadata_osaka.csvによるとサンプル１－１８までがhigashihori、19－32までがhonmachiのサンプルです。higashihoriとhonmachiの間で炭素分解能力に違いはあるでしょうか？この問いに答える方法を学ぶのがこのセクションの目標です。具体的には、１）その違いを一瞬で理解できるように可視化する方法と、２）その違いがランダムな過程から生まれた偶然かどうかを検定する方法、の二つを解説していきます。

いきなり31次元のデータを使うと、直感的イメージを得ることができないので最初の2列V1,V2だけを使って可視化の方法を考えます。まずは上の図と同様、棒グラフにしてみます。

次は、Rを使って散布図を描いてみましょう。基本はplot(x,y)でｘ－ｙの散布図が描けます。plot()関数にはいろいろオプションがあって図のデザインを変更できますが、ここでは各点の色をmetadata_osaka$place_sampleの値("higashihori"か"honmachi")によって、区別し、cexというオプションの値で各点の大きさを指定しています。

as.numeric(metadata_osaka$place_sample)

plot(osaka_summary_integ_ave$V1,osaka_summary_integ_ave$V2,col=2*as.numeric(metadata_osaka$place_sample),cex=2)

これの実行結果は下のようになります(見やすいようにWidthとHeightを変えています)。このグラフから分かることは、１）V1とV2の値はこの二次元プロット上にランダムに分布しているというよりも、左下から右上に向かっていくパターンで分布していること、さらに２）higashihoriからのサンプル（赤い丸）とhonmachiからのサンプル（青い三角）が混ざっていないことです。この二つの特徴はグラフを見れば一目瞭然ですね。これが二次元散布図のパワーです。すごい！

しかし、思い出してみると元のデータの次元は31です。31次元散布図なんて描けませんし、天才数学者か物理学者じゃなきゃ頭の中でも想像できません。なんとか31次元のデータを二次元散布図に落とし込めたら・・・。この欲求を満たす方法を、統計学や機械学習では「（教師なし）次元削減」といいます。

人間は2次元散布図は直感的に理解できるのでこれ以上シンプルにする必要はないのですが、まずはこの方法のアイデアの核を上の二次元散布図を使って説明します。この図の分布の傾向として、前述したように左下から右上に向かったパターンが見れるおかげで、次のような軸を設定すれば二次元を一次元に削減できます。

心眼をつかえば、なんとなく以下のように軸（黒い太線矢印）が見えてくるはずです。ここで注意していただきたいのは、この軸の方向は普通の回帰直線の方向とは一致しません。図にある点線のように各点から垂直に軸に線を引いたとき、すべての点についてこの線の長さの総和（＝軸と各点との距離の総和）が最小になるように方向を決めます。これ以上詳しい説明は統計の教科書やウェブに載っているので省略します。

この黒い軸上に点を移動させる（＝写像を取る）と、以下のような１次元の図になります。これが２次元から１次元への次元削減で、この軸を「主成分(Principal Component)」といい、このような解析を主成分分析(PCA: Principal Component Analysis)と呼びます。後ほど改めてスクリプトは説明するのでここでは省略しますが、以下のスクリプトを実行するとこのグラフが描けます。

library(vegan)

osaka_test0.d<-vegdist(osaka_summary_integ_ave[,1:2], method="euclidian")

PCoA_osaka0 <- cmdscale(osaka_test0.d, k = 1)

x<-PCoA_osaka0[,1]

y<-rep(0:0, length=32)

plot(x,y, cex=2, col=2*as.numeric(metadata_osaka$place_sample),pch=as.numeric(metadata_osaka$place_sample))

問題は、２次元の図から1次元の図にして、ビジュアルが簡単になった分、実は情報が失われていることです。上の２次元散布図で水色で囲った赤丸は３つあったのですが、それは主成分を取った１次元の散布図では一番左の三角のすぐ左の二つの点になってしまっています。2次元ではたしかに区別された二つの点はこの主成分上では重なってしまうわけで、その分情報が失われているということになります。そのため、PCAは、次元を削減したときに元の何パーセントの情報が残されているかという計算ができるのですが、ここではこれ以上説明しません。

さて、2次元から1次元への次元削減だけではイメージがつかみにくいと思うので、今度はosaka_summary_integ_aveのV1,V2,V3の3つのデータを使ってPCAの威力を確認したいと思います。単純に3次元散布図は以下のようなスクリプトで描けます。今度も3次元空間全体にサンプルが分布しているわけではないのが見て取れると思います。

#3D scatter plot

library(scatterplot3d)

scatterplot3d(osaka_summary_integ_ave$V1,osaka_summary_integ_ave$V2,osaka_summary_integ_ave$V3,color=2*as.numeric(metadata_osaka$place_sample),pch=1, cex.symbols=2,angle=50)

ここで再び心眼を使えば、以下のような水色の平面状に点が落としこめる気がしてきます。

心眼でやっても仕方ないのでちゃんとPCAのスクリプトを書いて（説明は省略）3次元データを二次元データに削減、プロットしてみましょう。グラフは下のようになります。

osaka_test.d<-vegdist(osaka_summary_integ_ave[,1:3], method="euclidian")

PCoA_osaka <- cmdscale(osaka_test.d, k = 2)

x<-PCoA_osaka[,1]

y<-PCoA_osaka[,2]

plot(x,y, cex=2, col=2*as.numeric(metadata_osaka$place_sample))

以上が主成分分析：PCAのイメージです。

Ecological Dissimilarityの計算：「似ている／似ていない・近い／遠い」を定量化する

さて、どうして２次元の散布図にすると、わかりやすいのでしょうか？ なぜなら経験上、近くにあるもの同士は似ている、遠くにあるもの同士は似ていないという感覚があるからでしょう。では、「似ている／似ていない」、「近い／遠い」ってなんでしょう？ この当たり前のことについてもう少し深く考えないとエコプレートを初めとして生態学でよく出てくる多変量の解析はできません。ここでは単純化した例を考えて見ます。サンプルが4つ（A,B,C,D）ありサンプルは2種類の基質V1,V2の分解能で特徴付けられているとします。A(3.0, 3.5), B(0.0, 0.5), C(1.0, 0.0), D(6.0, 6.0) という仮想的な値を与えます。

################Ecological Distance

A<-c(3.0,3.5)

B<-c(0.0,0.5)

C<-c(1.0,0.0)

D<-c(6.0,6.0)

test_sample<-as.data.frame(rbind(rbind(rbind(A,B),C),D))

colnames(test_sample)<-c("V1","V2")

View(test_sample)

単純に2次元散布図を描くスクリプトは以下のようになります。

#Simple 2d scatterplot

plot(test_sample, cex=0.5,xlim=c(0,7),ylim=c(0,7))

text(test_sample$V1, test_sample$V2, labels=rownames(test_sample), cex=0)

グラフは以下のようになります。

上図から人目でわかるように、BとCの距離が一番短く、それは分解能が、一番「似ている」という解釈になると思います。具体的に各サンプル間の距離を計算するには、veganパッケージ内のvegdist()という関数を使います。一般に日常会話も含めて「距離」と呼んでいるものは、「ユークリッド距離(Euclidean Distance)」であり、オプションでmethod="euclidean")を指定することで以下のように計算できます。

#examples of ecological dissimilarity

vegdist(test_sample, method="euclidean") #Euclidean

その結果は以下のようになります。

A B C

B 4.242641

C 4.031129 1.118034

D 3.905125 8.139410 7.810250

上の計算は、「ユークリッド距離によって、サンプル間の分解能の類似度を評価」する計算である、と言えます。あるいは、「（ユークリッド）距離が大きいほど類似度は小さい」、または「（ユークリッド）距離が大きいほど非類似度が大きい」とも言えます。ここで問題は、この計算による類似度の評価が、われわれの生物学もしくは化学の理解の観点から妥当かどうか、という点です。A(3.0, 3.5), B(0.0, 0.5), C(1.0, 0.0), D(6.0, 6.0) という値をひとつひとつ吟味してみると、AとDは基質V1とV2を共に分解できているがその強度・規模(magnitude)が異なるのに対し、BはV1のみ、CはV2のみしか分解できていません。このとき、妥当な解釈は、AとDの類似度が高く（非類似度が小さく）、BとCの類似度が小さい（非類似度が大きい）というものになるはずです。しかしユークリッド距離は以下のような定義で計算するので、AとDの各要素（V1,V2)のそれぞれの絶対値が大きく反映され、非類似度が大きくなってしまいます。

そこで数多くの「非ユークリッド」距離が定義されています（全貌は、Tree diversity analysis manualの第8章参照）。その多くは計量性 (metric)の条件の一部を満たしていないのですが（「計量性」説明はwikipedia参照）、より上で述べた解釈をより強く反映できます。

ここでは、まず上のA-Dのデータを値がゼロより大きいか多くないかの２値に変換後（すなわち、存不在:presence/absence）に非類似度を計算する、Jaccard距離を以下のスクリプトで計算してみます。技術的な注意点ですが、jaccard距離は一般には連続値データに対しても定義されるのでオプションでbinary=TRUEを指定しないと意図しない値が出てしまいます。

#Non-Euclidean (& non-metric or semimetric) dissimilarity

vegdist(test_sample, method="jaccard", binary=TRUE) #Presence/absence

この結果は以下のようになります。

A B C

B 0.5

C 0.5 1.0

D 0.0 0.5 0.5

見て分かるとおり、AとDのjaccard距離が一番小さくゼロ、共通の分解能を持たないBとCの間の距離が1.0となって最大になりました。数学的な注意点としては、AとDはもともと同一のデータではないにもかかわらず、距離がゼロと計算されるため、jaccard距離は計量性の一部を満たしていません（ユークリッド距離は計量性を満たしており、距離がゼロとなるのは、同一データのときのみです）。参考までに二値の場合のjaccard距離の定義を以下に示します。

もうひとつ、Jaccard距離と似た挙動をしめしつつデータの連続値も加味したものとして、Bray-Curtis距離を計算してみます。Jaccard距離と同様に、B,C間の距離が最大、A,D間の距離が最小となっている一方、Jaccard距離との違いは例えば、A,B間とD,B間の距離が異なっています。

> vegdist(test_sample, method="bray") #Bray-curtis

A B C

B 0.8571429

C 0.7333333 1.0000000

D 0.2972973 0.9200000 0.8461538

参考までにBray-Curtis距離の定義を以下に示します。

主座標分析(PCoA: Principal Coordinate Analysis)の基礎

上で紹介した非ユークリッド距離を用いて非類似度を計算しそれを可視化したい場合は、主成分分析は使いません。なぜなら、主成分分析では情報の損失量を計算する過程でユークリッド距離が定義される世界でしか成り立たない三平方の定理を用いているからです。主成分分析の代わりになるのが主座標分析と呼ばれる手法で英文略称でPCoAとも言います。計算方法の詳細は説明しませんが、PCoAでは、PCAのように元のデータから出発するのではなく、各サンプル間の距離の情報（距離行列：vegdistで計算したもの）を使って、サンプル間の距離の値を維持したまま、（ユークリッド距離が定義されている）ユークリッド空間上に各サンプルを配置しなおすという計算をします。

技術的な注意

PCoAは、非ユークリッド距離で定義された各サンプル間の距離をユークリッド距離の世界に写しこむ作業なので当然ゆがみが生じる点は注意が必要です。つまり、通常は非常に多次元なデータを二次元の散布図にするため、二次元の散布図上での各サンプル間の距離（ユークリッド距離）は元の非ユークリッド距離に基づく距離の値とは一致しません。ただし、ｎ次元のデータをｎ次元のユークリッド空間に配置した場合には、非ユークリッド距離で定義された距離の値は、ｎ次元ユークリッド空間上のユークリッド距離としてゆがむことなく完全に再現されます。この意味でPCoAは距離を維持した手法のため、計量的多次元尺度法(metric Multi-Dimensional Scaling: metric MDS）とも呼ばれます。計量性の呪縛を断ち切り、距離の違いを非類似度のランクに変換してから可視化する方法は非計量的多次元尺度法、non-metric MDS, NMDSとよばれていますが、ここでは時間の関係上説明しません。以下のウェブログの解説が詳しいです。

https://jonlefcheck.net/2012/10/24/nmds-tutorial-in-r/

さて、実際にA,B,C,Dのデータを使って可視化してみましょう。

まずは在不在に基づく非類似度であるJaccard距離によるプロットです。これには、（１）vegdist()関数で距離行列を計算し、（２）cmdscale()にその距離行列を引数として指定、二次元の可視化を行うとしてオプションをk = 2に指定してPCoAを実行します。その後、（３）PCoAの結果から各点の二次元ユークリッド空間上での座標値(x, y)を取得し、plot()とtext()を組み合わせて二次元散布図を描きます。スクリプトは以下の通りです。

#PCoA based on Jaccard distance

test_ja.d<-vegdist(test_sample, method="jaccard", binary=TRUE) #calculate Jaccard dissimilarity

PCoA_test_ja <- cmdscale(test_ja.d, k = 2) #execute PCoA

x<-PCoA_test_ja[,1] #extract the position of each sample (A,B,C,D)

y<-PCoA_test_ja[,2]

plot(x,y, cex=0.5,xlim=c(-0.5,0.5),ylim=c(0,0.1)) #2D scatter plot

text(x,y,label=rownames(test_sample)) #set the symbols

可視化の結果は以下の通りです。Y軸方向には値の差が小さく、ほぼ一直線に並んだ結果です。さらにAとDの距離はゼロなので重なっています。

次にBray-Curtis距離を用いたプロットです。上のスクリプトを少しだけ変更すればよく、以下のようになります。

#PCoA based on Bray-Curtis distance

test_bc.d<-vegdist(test_sample, method="bray") #calculate Bray-Curtis dissimilarity

PCoA_test_bc <- cmdscale(test_bc.d, k = 2) #execute PCoA

x<-PCoA_test_bc[,1] #extract the position of each sample (A,B,C,D)

y<-PCoA_test_bc[,2]

plot(x,y, cex=0.5,xlim=c(-0.7,0.4),ylim=c(-0.7,0.4)) #2D scatter plot

text(x,y,label=rownames(test_sample)) #set the symbols

図は以下の通りです。Jaccard距離に基づく結果とはだいぶ違います。ただ、基本的に距離の大小の関係は似ています。

最後に、ユークリッド距離に基づくPCoAの図です。図の縦横比（アスペクト比）を１：１にすれば、もともとの二次元散布図とほぼ同じ配置になっているのがわかるはずです。

#PCoA based on Euclidean distance

test_eu.d<-vegdist(test_sample, method="euclidean") #calculate Euclidean dissimilarity

PCoA_test_eu <- cmdscale(test_eu.d, k = 2) #execute PCoA

x<-PCoA_test_eu[,1] #extract the position of each sample (A,B,C,D)

y<-PCoA_test_eu[,2]

plot(x,y, cex=0.5,xlim=c(-3,6),ylim=c(-3,6)) #2D scatter plot

text(x,y,label=rownames(test_sample)) #set the symbols

さて上のユークリッド距離に基づくPCoA散布図ですが、実はこれは主成分分析PCAの結果と全く同じになります。Rのveganパッケージに含まれる関数rda()でPCAを行うには通常以下のようなスクリプトを用意します。scaling=1というオプションでの可視化では、サンプル間の距離は、サンプル間のユークリッド距離に対応します。

##########Finally, show PCA with scale 1 (see Biodiversity manual Chapter 10)

PCA_test<-rda(test_sample)

plot1<-ordiplot(PCA_test, scaling=1, type="text")

ｘ軸の向きが上のPCoAと逆転していますが、各点の相対的位置はほぼ同一であることが読み取れると思います。

PC軸が何を示すとか、説明率とか詳細はいろいろあるのですが、ここでは触れずに終了とします。