Random Forest

十〜数百、それ以上のデータがある状態で分類器を構築する場合、私はRandom Forestをよく使います。Random Forestとは、標本データを復元ありの無作為抽出 (bootstrap)して作成した仮想データを多数生成して、それぞれのデータに対して毎回ランダムに選択した変数群を用いて決定木を構築、各々の決定木の多数決で予測を行うといった分類器です。bootstrapでデータの揺らぎを学習し (bagging)、多数の異なる決定木の多数決 (ensemble)でモデルの揺らぎを学習するイメージです。予測精度が高く、過学習 (overfit)しない点、bootstrapにより予測率が評価されるので、cross validation等が必要ない点など、扱いやすい方法です。 Rで計算できます。

パッケージをロードして、データを読み込みます (データはここからDLできます)。
library('randomForest')
tr<-read.table("./lung_michigan.txt",header=T,row.names=1)

Random Forestを実行します。パラメータは２つ。ntreeは何個の決定木を構築するか (=bootstrap数) 、mtryは決定木のサイズです 。※mtryはデフォルトでは変数の数の平方根です。sampsizeは省略できますが、群のサイズがアンバランスな場合、多い群に予測が引っ張られるため、小さい群のサイズを指定します。
rf<-randomForest(class~.,data=tr,importance=T,proximity=T,ntree=1000,mtry=10,sampsize=c(24,24))

このデータはA:alive=24, D:death=62で構成されています。予測精度は以下で確認できます。
rf$confusion
A D class.error
A 57 5 0.08064516
D 9 15 0.37500000

ntreeに対する予測精度の推移をプロットできます。
plot(rf)


緑が全体の予測精度（error rate）、赤がA、黒がDの予測精度です。

予測に対する各変数の寄与度を得る事ができます。
imp<-importance(rf,scale=TRUE)
imp <- imp[, -(1:(ncol(imp) - 2))]
imp
MeanDecreaseAccuracy MeanDecreaseGini
SEC31L1 1.0843264667 0.90073387
RAFTLIN 0.6662848410 0.54427994
SLC2A1 0.4065791941 0.21934335
・・・・

寄与度のプロットも可能です。
varImpPlot(rf)


寄与度の計算は決定木を構築する際、該当変数をモデルから除いた際の、予測精度の低下 (Mean Decrease Accuracy)、あるいはGini indexの減少 (Mean Decrease Gini)に基づいています。

Random Forestでは、予測モデルの構築と同時にサンプル間の類似性も計算することが可能です。距離行列に変換してWard法でクラスタリングしてみます。
ds<-as.dist(1.0-rf$proximity)
hc<-hclust(ds,"ward")
plot(hc)

多重補正 (その３)

FDRの推定法によく用いられる"BH"などはtail area-based FDR (FDRと表記) と呼ぶようです。FDRはp値がα未満のときのfalse positiveの割合をあらわします。一方、Efronらがlocal FDR (fdrと表記) を提案していて、これはp値がαのときのfalse positiveの確率をあらわすそうです。ともあれ、関心は有意と考えている仮説群の確からしさで、どのくらい嘘が混じるリスクがあるかです。Rのfdrtoolパッケージでは、統計量のリストからempiricalにfalse positiveの割合etaを計算してくれます。以下は、500のp値リストを用いたシミュレーションです。数千〜数万の検定が発生するゲノムスケールの解析に有用かと。

library("fdrtool")
x<-rnorm(500,m=c(rep(0,250),rep(3,250)))
p<-2*pnorm(sort(-abs(x)))
fdrout<-fdrtool(p,statistic="pvalue")



推定されたパラメータは
fdrout$param

p値のリストは
fdrout$pval

local FDRは
fdrout$lfdr

tail area-based FDRは
fdrout$qval

で参照できます。

しょぼEXCEL

BMC Bioinformatics. 2004 Jun 23;5:80.
Mistaken identifiers: gene name errors can be introduced inadvertently when using Excel in bioinformatics.

Excelで特定の遺伝子名を含むファイルを読み込むと該当遺伝子名が自動的に変換されてしまうという不具合に関するレポート。例えば、以下のような遺伝子名が日付・実数に変換されてしまう。

DEC-1 -> 1-Dec
2310009E13 -> 2.31E+13.

後者はRiken ID。遺伝子名でIDマッチしている類いのシステム (GSEA MSigDBとか)を使うときは注意。しかし、余計なことを...。

Now available: MSigDB v2.5

GSEAで用いる遺伝子セットのデータベースMSigDBのv2.5がリリースしたようです。最新のMSigDBはC1〜C5のカテゴリで公開されています。

C1：染色体上の位置
C2：パスウェイ (新たにKEGGが追加)
C3：転写モチーフ
C4：共発現 (新たにがんにおける共発現モジュールが追加)
C5：Gene Ontlogy (新設)

GSEA初期の頃に比べ、大幅にコンテンツが増えました。

多重補正 (その２)

古典的なアプローチによる多重補正の計算は比較的簡単に計算できます。プログラミングも容易です。統計パッケージRにももちろん多重補正のパッケージが存在します。p.adjustを使って多重補正をシミュレーションしてみました。
x<-rnorm(50,m=c(rep(0,25),rep(3,25)))
p<-2*pnorm(sort(-abs(x)))
holm<-p.adjust(p,"holm")
hommel<-p.adjust(p,"hommel")
hochberg<-p.adjust(p,"hochberg")
bonferroni<-p.adjust(p,"bonferroni")
BY<-p.adjust(p,"BY")
BH<-p.adjust(p,"BH")
fdr<-p.adjust(p,"fdr")
pmat<-cbind(p,bonferroni,hochberg,hommel,BY,BH,fdr)
matplot(pmat,pch=1:ncol(pmat),type="o",col=rainbow(ncol(pmat)),main="comparison of multiple test correction",xlab="position of p-value in ascending order",ylab="corrected p-value")
legend(1,max(pmat),legend=colnames(pmat),col=rainbow(ncol(pmat)),pch=1:ncol(pmat))



ランダムに生成した50回の検定によるp値 (赤丸でプロット) に対して、各種補正法を適用しプロットしたものです。"BY", "BH"はFWERではなくFDR (False Discovery Rate) をコントロールする手法です。FDRはFWERより検出力の高い手法となります。FDRの計算は、変数が非常に多いケースにおいて、コンピュータの速度向上も伴いempiricalな方法も登場してきています。この辺はまたの機会で。

Hochberg, Y. (1988). A sharper Bonferroni procedure for multiple tests of significance. Biometrika, 75, 800–803.
Hommel, G. (1988). A stagewise rejective multiple test procedure based on a modified Bonferroni test. Biometrika, 75, 383–386.
Benjamini, Y., and Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. Annals of Statistics 29, 1165–1188.
Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B, 57, 289–300.