EcoPlate解析 (４)

６）多変量解析の基礎２

ここでは、エコプレートデータを用いて、PCoAによる可視化・および処理間の違いの検定を行います。

EcoPlate解析（３）ではEcological Dissimilarityを簡単に説明しましたが、異なるDissimilarityがどのような意味付けになるかは、Anderson et al. (2011) Navigating the multiple meanings of beta diversity: a roadmap for the practicing ecologist. Ecology Letters 14:19-28の図５を確認する必要があります。

まずは、今後の準備として、エコプレートに含まれる31種類の炭素基質の名前を読み込みます。

#Prepare the list of consistance name of chemicals included in ecoplate

chem_name<-c("Pyruvic-Acid-Methyl-Ester", "Tween-40", "Tween-80","alpha-Cyclodextrin", "Glycogen", "D-Cellobiose","alpha-D-Lactose", "beta-Methyl-D-Glucoside", "D-Xylose", "i-Erythritol", "D-Mannitol","N-Acetyl-D-Glucosamine", "D-Glucosaminic-Acid", "Glucose-1-Phosphate","alpha-Glycerol-Phosphate","D-Galactonic-Acid-gamma-Lactone", "D-Galacturonic-Acid", "2-Hydroxy-Benzoic-Acid", "4-Hydroxy-Benzoic-Acid", "gamma-Hydroxybutyric-Acid", "Itaconic-Acid", "alpha-Ketobutyric-Acid", "D-Malic-Acid", "L-Arginine", "L-Asparagine", "L-Phenylalanine", "L-Serine", "L-Threonine", "Glycyl-L-Glutamic-Acid", "Phenylethyl-amine", "Putrescine")

それから、osaka_summary_integ_ave (時間積分し繰り返しの最大値を採ったもの)を例に、行と列にそれぞれ名前をつけます。

rownames(osaka_summary_integ_ave)<-metadata_osaka$sample_ID

colnames(osaka_summary_integ_ave)<-chem_name

次に、前処理としてマイナスの値についてゼロへと変換します。マイナスの値はウェルの発色がコントロールよりも低かったことを意味するのでゼロへの変換し、対応する基質に対する分解能力が全く無かった、と解釈することにします。

#first convert negative values into zero (negative values imply that the color development was weaker than control)

osaka_summary_integ_ave[osaka_summary_integ_ave<0]<-0

PCoAによるエコプレートデータの可視化

つづいて、Euclidean distanceとBray-Curtis distanceをそれぞれ計算します。

#calculate Euclidean distance

osaka_test_eu.d<-vegdist(osaka_summary_integ_ave, method="euclidean")

#calculate Bray-Curtis distance

osaka_test_bc.d<-vegdist(osaka_summary_integ_ave, method="bray")

Bray-Curtis distanceを計算するスクリプトを実行すると、警告メッセージがコンソールに表示されます。これはすべての値がゼロとなるサンプル(行)が含まれるのが原因です。

> osaka_test_bc.d<-vegdist(osaka_summary_integ_ave, method="bray")

Warning message:

In vegdist(osaka_summary_integ_ave, method = "bray") :

you have empty rows: their dissimilarities may be meaningless in method “bray”

そこでその行をデータおよびメタデータから削除します（2行目のサンプルです）。

#Exclude samples with all zero values from data and metadata

osaka_summary_integ_ave<-osaka_summary_integ_ave[-2,]

metadata_osaka<-metadata_osaka[-2,]

その上でもう一度、Bray-Curtis distanceを計算します。今度は警告が出ないはずです。

#calculate Bray-Curtis distance

osaka_test_bc.d<-vegdist(osaka_summary_integ_ave, method="bray")

こんどはちょっと別の距離も計算してみます。EcoPlate解析（３）の冒頭のエクセルのグラフを見ればわかるように、サンプル１３以降はそもそも全体的に値が大きいことが分かります。この値の絶対値の影響を簡単に取り除くには、生データを割合データに変換すればよいです。veganパッケージにはdecostand()という関数が用意されていて、割合などへの変換によって生データを標準化することができます。割合を計算するにはmethod="total"を指定します。

#Calculate proportion

osaka_summary_integ_ave_prop<-decostand(osaka_summary_integ_ave, method="total")

この標準化したデータでEuclidean distanceを計算してみます。

#calculate Euclidean distance on proportion

osaka_test_prop_eu.d<-vegdist(osaka_summary_integ_ave_prop, method="euclidean")

これらの距離それぞれを使ってPCoAを実行するには以下のようなコードを書き、順々にグラフへと可視化できます。

#PCoA for each distance

PCoA_osaka_eu <- cmdscale(osaka_test_eu.d, k = 2)

x<-PCoA_osaka_eu[,1]

y<-PCoA_osaka_eu[,2]

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

PCoA_osaka_prop_eu <- cmdscale(osaka_test_prop_eu.d, k = 2)

x<-PCoA_osaka_prop_eu[,1]

y<-PCoA_osaka_prop_eu[,2]

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

PCoA_osaka_bc <- cmdscale(osaka_test_bc.d, k = 2)

x<-PCoA_osaka_bc[,1]

y<-PCoA_osaka_bc[,2]

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

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

この可視化結果を見ると、だいぶ距離の定義によってサンプルの分布ぐあいが違います。まずEuclidian distanceを使った二つの結果とBray-Curtis distanceを使った結果の違いの原因は二つあります（冒頭で紹介した論文参照）。まず、Euclidean distanceでは2サンプルで共に値がゼロとなる項目(列)も類似度を高めることになる一方(joint absences are considered)、Bray-Curtis distanceではそのような項目は無視されます。次に、EcoPlate解析（３）でも確かめたように、Euclidean distanceが生データの絶対値を強く反映するのに対し、Bray-Curtis distanceはより組成(composition)と相対頻度(relative abundance)を反映します。Euclidean distanceが絶対値を強く反映する性質を弱めたものが、proportionを最初に計算してからEuclidean distanceを計算したものになります。したがって、Euclidean on proportionの結果はよりBray-Curtis on raw dataに近い結果になっています。

エコプレートパターンに処理間で違いはあるか (db-RDA)？

つぎに、検定の話をします。上のPCoAのグラフを見れば、サンプリング場所の異なる丸(higashihori)と三角(honmachi)の間になんとなく差があるように見えます。しかし丸印のサンプル間にも、三角印のサンプル間にもばらつきがあるので、このようなばらつきが本当は場所間に差が無いのに偶然生じたのかどうか判断できません。また、データ間の全体のばらつき（＝丸印間のばらつき＋三角印間のばらつき＋丸－三角間のばらつき）のうち、何%程度が場所間（あるいは一般には処理間）のばらつきなのか分かりません。あるいは、場所のようなカテゴリー変数ではなく、たとえば温度のような連続値変数の違いによってばらつきがどの程度説明できるか知りたい、という状況はよくあります。

説明したい変数がエコプレートのデータのように多変量ではなく、たとえば種数や一次生産といった単変量の場合には、分散分析 (ANOVA)や回帰分析を行えばよいわけです。これが多変量がターゲットになった場合は、主成分分析(PCA)もしくは主座標分析(PCoA)と回帰分析を組み合わせたRedundancy Analysis (RDA)もしくはdistance-based Redundancy Analysis (db-RDA)という手法を用います。カテゴリー変数間での差をみるにはは、PERMANOVA(NPMANOVA)およびPERDISPという方法を使うことが多いのですが、簡単のためここでは説明せず、distance-based RDAを用います。詳細はリンク先を参照ください。また、veganパッケージのFAQも大変参考になります。

db-RDAはPCAではなくPCoAを基にしているのでユークリッド距離以外にも使えます。db-RDAを実行するには、capscale()という関数を使います。関数の()の中には、「統計モデル式」を書く必要があります。まずはユークリッド距離で処理したデータを使いますが、ここでの目標は、サンプル間の炭素分解能力の違い(＝ユークリッド距離: osaka_test_eu.d)がサンプル場所間(metadata_osaka$place_sample)で異なるかどうか、を知ることですので、統計モデル式は、 "目的変数(response variable)～説明変数(explanatory variable)"という形式で以下のように書きます（"従属変数(dependent variable)～独立変数(independent variable)"という言い方もあります）。

osaka_test_eu.d~metadata_osaka$place_sample

実際にdb-RDAを実行しその結果を新しいオブジェクトosaka_test_eu_rdaに格納するには以下のスクリプトを実行します。

##distance-based RDA analysis

osaka_test_eu_rda<-capscale(osaka_test_eu.d~metadata_osaka$place_sample)

この結果osaka_test_eu_edaの一部を使って説明変数で説明されるばらつきの割合は以下のように計算できます。計算された値は0.4530333です。

#The variance explained (constrained) by explanatory variable(s)

osaka_test_eu_rda$CCA$tot.chi/(osaka_test_eu_rda$CCA$tot.chi+osaka_test_eu_rda$CA$tot.chi)

計算の方法が違うので厳密には同じ値になりませんが、db-RDAに含まれるPCoAの計算で出てくる固有値を使っても、この説明割合を計算できます。

#Another way of calculation

summary(eigenvals(osaka_test_eu_rda))

この結果はコンソールに固有値の一覧として表示されますが、説明変数（ここではサンプル場所）によって説明される割合は、Cumulative ProportionのCAP1の値です。一般に説明変数がn個ある場合は、説明変数によって拘束される軸がn個出てくるので、CAPnの真下のCumulative Proportionの値を見ればよいことになります。

> summary(eigenvals(osaka_test_eu_rda))

Importance of components:

CAP1 MDS1 MDS2

Eigenvalue 0.2696 0.2208 0.02764

Proportion Explained 0.4530 0.3710 0.04644

Cumulative Proportion 0.4530 0.8240 0.87049

Permutationによる検定

さて、上の結果では約45%のばらつきは、サンプル場所の違いにより説明できることがわかりました。さてこの値は「サンプル場所間でエコプレートパターンに差は無い」という帰無モデルからは生じ得ないのでしょうか？ 差は無いと仮定したモデルからも45%とか55％とかの値が出る場合には、45%という値が出たからといって「差が無い」という帰無モデルは棄却できません。単変量の解析である分散分析では、自由度あたりの説明変数によるバラつきと自由度あたりのその他のバラつきの比を取ってF値 (F value) と呼び、その値が大きいほど説明変数による効果が大きいことを示します。古典的な分散分析では母集団が正規分布を持つと仮定してF値の確率分布を数学的に計算できるため、F値の実測値とその確率分布を利用して、「差が無い」という帰無モデルからF値の実測値よりも大きな値が出る確率Pを計算できます。この確率Pが十分低い (たとえばP<0.01) とき、帰無モデルからデータが生じた可能性は低いとして帰無モデルを棄却することになります。これがいわゆる検定です。たとえば、上のdb-RDAの結果についても擬似的なF値 (pseudo-F value)を計算可能で、以下のスクリプトによると値は約24.01となります。

#Check pseudo-F value (see F value in ANOVA)

df_re<-length(osaka_summary_integ_ave) - 1 - 1 #degree of freedom

osaka_test_eu_rda$CCA$tot.chi/(osaka_test_eu_rda$CA$tot.chi/df_re)

ここで問題は、一般に生態学で扱うような多変量は多次元正規分布には従わないことです。したがってpseudo-F valueが従う確率分布は既知ではありません。このような状況下では一般的に"permutation"という方法で手元にあるデータからランダムにデータを再生成し検定を行います。正確さを犠牲にして具体例で簡単に説明すると、「サンプル場所間で違いが無い」という帰無モデルのデータを生成するには、もともとのデータosaka_summary_integ_aveのサンプル場所ラベルmetadata_osaka$place_sampleをシャッフルすればよいわけです。ラベルのシャッフルは以下のようにsample()という関数を使って行うことができます。

####Example of permutation (shuffling)

set.seed(123) #fix the random seed

metadata_osaka_perm<-metadata_osaka

leng<-length(metadata_osaka[,1])

z<-sample(1:leng, leng, replace=FALSE) #re-sampling without replacement

z

for(i in 1:leng) metadata_osaka_perm$place_sample[i]<-metadata_osaka$place_sample[z[i]]

最後のfor loopでラベルの入れ替えを行っています。この入れ替えたラベルに基づいてEuclidean distanceを用いてPCoAとdb-RDAを実行しF値を計算するには以下のスクリプトを実行します。

#plot permutated data by PCoA

PCoA_osaka_eu <- cmdscale(osaka_test_eu.d, k = 2)

x<-PCoA_osaka_eu[,1]

y<-PCoA_osaka_eu[,2]

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

#db-RDA for permutated data

osaka_test_eu_perm_rda<-capscale(osaka_test_eu.d~metadata_osaka_perm$place_sample)

#Check pseudo-F value (see F value in ANOVA)

osaka_test_eu_perm_rda$CCA$tot.chi/(osaka_test_eu_perm_rda$CA$tot.chi/df_re)

コンソールに表示されるF値の値を記録した後、z<-sample(1:leng, leng, replace=FALSE) #re-sampling without replacement以降のスクリプトを繰り返してもう一度シャッフルしF値を計算してみます。この結果は以下のグラフにまとめています。ラベルが置き換わっているのと、グラフに重ねてF値が計算されています。Permutation後のデータのF値は実測値（original data）よりもずっと小さいのがわかります。

上の説明はpermutationの概念を理解するため詳細を省いているのに注意していただきたいのですが、これで概念が理解できたら、anova()という関数を使ってpermutationを好きな数だけ繰り返し、permutationからのF値がオリジナルのF値を越える回数を計算することができます。スクリプトはとても単純で以下のようになります。その結果はコンソールに表示されます。

#permutation test

anova(osaka_test_eu_rda, permutations=1999, by='terms')

その結果はコンソールに表示されます。F値がオリジナルを超える確率 Pr(>F)が5e-4と非常に低いのが分かるでしょう。

> anova(osaka_test_eu_rda, permutations=1999, by='terms')

Permutation test for capscale under reduced model

Terms added sequentially (first to last)

Permutation: free

Number of permutations: 1999

Model: capscale(formula = osaka_test_eu.d ~ metadata_osaka$place_sample)

Df SumOfSqs F Pr(>F)

metadata_osaka$place_sample 1 0.26962 24.02 5e-04 ***

Residual 29 0.32553

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Euclidean distanceについて計算したのと同様、割合に対するEuclidean distanceと、Bray-Curtis distanceを用いた結果についてもdb-RDAとpermutation testを以下のように実行してみます。

#db-RDA analysis

osaka_test_prop_eu_rda<-capscale(osaka_test_prop_eu.d~metadata_osaka$place_sample)

osaka_test_prop_eu_rda$CCA$tot.chi/(osaka_test_prop_eu_rda$CCA$tot.chi+osaka_test_prop_eu_rda$CA$tot.chi)

anova(osaka_test_prop_eu_rda, permutations=1999, by='terms')

osaka_test_bc_rda<-capscale(osaka_test_bc.d~metadata_osaka$place_sample)

osaka_test_bc_rda$CCA$tot.chi/(osaka_test_bc_rda$CCA$tot.chi+osaka_test_bc_rda$CA$tot.chi)

anova(osaka_test_bc_rda, permutations=1999, by='terms')

その結果、Euclidean distance on proportion, Bray-Curtis distanceを用いたときのサンプル場所間でのばらつきの説明率（=statistical power）は、それぞれ、12.6%, 29.1% であることが分かります。以上がエコプレートを例にとった多変量解析の例です。

炭素分解「多機能性」に処理間で違いはあるか (quantile-based multifunctionality)？

この最後のセクションではエコプレートパターンという多変量をmultifunctionalityという単変量に変換する方法について解説します。multifunctionalityの計算手法には複数あるのですが、詳しくはわれわれの論文(T. Miki, T. Yokokawa, K. Matsui* (2014) Proceedings of The Royal Society B vol.281 no.1776 20132498 )の中で引用されている文献を参照してください。

さて、quantile-based multifunctionalityという計算方法では、各機能（＝各基質の値）ごとに順位順に並べ、下から数えてX％については機能無し、それ以上は機能有りと二値化(binarization)する方法です。コアとなる関数はquantile()です。

たとえば、osaka_summary_integ_aveのデータを新しいデータフレームにコピーした後、4つ目の基質chem_name[4]に注目して、下から90％の値を計算してみるには以下のようなスクリプトを実行します。

#Convert data into binary values by quantile (quantile-based multifunctionality)

file_summary<-osaka_summary_integ_ave

i<-4

chem_name[4]

thres<-0.9

quantile(file_summary[,i], thres)

この方法をfor loopの中で使えばデータ全体を二値化できます。

mf<-list() #prepare the output list

thres<-0.9

for(i in 1:length(file_summary)) mf[[i]] <-(file_summary[,i] > quantile(file_summary[,i], thres)) # return true if the color value qunatile is greater than threshold

M_mf<-as.data.frame(mf) #converted into dataframe

row.names(M_mf)<-metadata_osaka$sample_ID

colnames(M_mf)<-chem_name

M_mf[M_mf==TRUE]<-1 #converting T,F to 1 and 0; ts is now dataframe to represent the binary matrix (E_BC in Figure 3)

View(M_mf)

この二値化されたデータについてapply()という関数を使って各行の合計値sumを計算すれば、サンプルごとのmultifunctionalityが0-31の間で計算できます。

#calculate multifunctionality

MF_osaka<-apply(M_mf,1,sum)

ハコひげ図（boxplot）を使ってサンプル場所間での比較を可視化するためにデータを分離した後、boxplot()を使ってみます。

#separate data for boxplot

higashihori_MF<-MF_osaka[metadata_osaka$place_sample=="higashihori"]

honmachi_MF<-MF_osaka[metadata_osaka$place_sample=="honmachi"]

boxplot(higashihori_MF, honmachi_MF, names=c("higashihori","honmachi"))

結果は以下のようになります。検定にかけたわけではありませんが、honmachiのほうがmultifunctionalityが高い傾向が見て取れます。これも多変量、とくに生態系機能に関するパラメータに関してよく使われる解析の一つです。