菌叢解析結果の比較 ( 3 ) 二元配置分散分析

前回の記事で、バートレット検定でphylumごとのrelative abundanceの分散が等しいかどうかを判定しました。p>0.05で分散が等しいと判定したphylumであるBacteria;__Actinobacteriaについて行います。データは引き続き前回のものを使います。Rで二元配置分散分析を行うときはaov関数を使用します。

[splus]
#.チルダ"~"の左側の変数を右側の式で説明するモデルを作成する
summary(aov(X[,3]~X$Test+X$Subject))
[/splus]

を実行すると

Df Sum Sq Mean Sq F value Pr(>F)
X$Test 4 0.01197 0.00299 0.826 0.54093
X$Subject 1 0.06030 0.06030 16.629 0.00277 **
Residuals 9 0.03263 0.00363

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

の様に返って来ます。つまり抽出法では有意差はありませんが、被検者間では有意差が出ました。これを(分散が等しいとして)すべてのphylumに行うべくスクリプトを書いてみます。

実は2要因で行うaov()には厄介な点があります。aov()をそのまま実行すると

[splus]
aov(X[,3]~X$Test+X$Subject)
[/splus]

Terms:
X$Test X$Subject Residuals
Sum of Squares 0.01197406 0.06029833 0.03263396
Deg. of Freedom 4 1 9

Residual standard error: 0.06021624
Estimated effects may be unbalanced

となってしまうので、summary(aov())の形で結果を抽出します。ここで欲しいのはP値ですから、Pr(>F)を呼び出したいのですがなかなかうまくいきませんでしたが、どうやら

[splus]
ASummary = summary(aov(X[,3]~X$Test+X$Subject))
ASummary[[1]]["Pr(>F)"][1,1] #.TestのP値
[1] 0.6532607
ASummary[[1]]["Pr(>F)"][2,1] #.SubjectのP値
[1] 0.5708131
[/splus]

とすれば、それぞれP値が返って来ます。それを踏まえると、

[splus]
ATable=NULL
for ( i in 3:ncol(X))
{
ASummary = summary(aov(X[,i]~X$Test+X$Subject))
AResult = c(colnames(X)[i],ASummary[[1]]["Pr(>F)"][1,1],ASummary[[1]]["Pr(>F)"][2,1])
ATable = rbind(ATable,AResult)
}
colnames(ATable)=c("Phylum","Test$P.value","Subject$P.value")
write.csv(ATable,"AOV_Result.csv",row.names=FALSE)
[/splus]

と書けば、少し回りくどいスクリプトですが、

“Phylum”,”Test$P.value”,”Subject$P.value”
“Bacteroidetes”,”0.540929440422521″,”0.00276660508206454″
“Firmicutes”,”0.703002431390229″,”0.333099501465002″
“TM7″,”0.687474809867674″,”0.44333185016966”

と、” AOV_Result.csv”にすべてのphylumに関して上手く結果が吐き出されます。有意差の判定は面倒臭いのでつけませんでした。ただ、実際は分散が等しい場合と等しくない場合で検定をかけなければいけないのですが、それはまた後の記事で書きます。

次の記事では、分散が等しくない場合に行うフリードマン検定について書きます。

Share on Facebook0Share on Google+0Share on Tumblr0Tweet about this on TwitterEmail this to someone

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です