2008-05-28
条件付き独立性の検定
連続値をとる変数間の (条件付き) 独立性の指標は (偏) 相関係数です (ただし多次元正規分布の仮定のもと)。偏相関係数の計算／検定法は割とよく知られているので、ここでは離散値 (カテゴリカル) をとる変数間の (条件付き) 独立性の検定法を紹介します。
帰無仮説が真のとき、以下のG2統計量はカイ2乗分布に従います。
ここで、N_ijkはセルijkの観測頻度、m_ijkはセルijkの期待頻度です。kはSのカテゴリーが取り得るすべてのパターンをなめます。自由度dfは(X_iのカテゴリー数-1)*(X_jのカテゴリー数-1)とSのすべての要素のカテゴリー数の積です。このG2統計量を用いて周辺・条件付独立性を検定するには、自由度dfのカイ二乗分布におけるパーセンタイルを考えればよいことになります。
帰無仮説が真のとき、以下のG2統計量はカイ2乗分布に従います。
ここで、N_ijkはセルijkの観測頻度、m_ijkはセルijkの期待頻度です。kはSのカテゴリーが取り得るすべてのパターンをなめます。自由度dfは(X_iのカテゴリー数-1)*(X_jのカテゴリー数-1)とSのすべての要素のカテゴリー数の積です。このG2統計量を用いて周辺・条件付独立性を検定するには、自由度dfのカイ二乗分布におけるパーセンタイルを考えればよいことになります。
2008-05-24
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)
パッケージをロードして、データを読み込みます (データはここから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)
2008-05-20
グラフライブラリ
グラフ構造に対するアルゴリズムを実装するには、グラフを表現するデータ構造を選択する必要があります。ノード間の関係は行列によりあらわすことができます。これは隣接行列と呼ばれ、例えば、i,j要素が1のときノードiとノードjの間にエッジがあり、0のとき無いといった具合に表現できます。当然、スパースなネットワークの場合、0だらけの行列となりメモリ効率が悪くなります。別の表現方法としてリスト構造があります。ポインタで数珠つなぎすることによって表現します。今は、グラフに関するライブラリが充実しているため、グラフ理論の多くのアルゴリズムの恩恵を受けることができます。私はC／C++系のコーディングではboost graph libraryを、perlではGraphモジュールをよく用います。例えば以下のような２項関係を記述したファイルがあるとします。
>example.txt
nodeA<tab>nodeB
nodeB<tab>nodeC
nodeC<tab>nodeA
この２項関係をグラフで表すと"nodeA->nodeB->nodeC->nodeA"とcyclicなグラフとなります。このようなファイルを入力としてDAG (Directed Acyclic Graph : 非巡回有向グラフ)かどうかを判定するプログラムは以下のように簡単に書けます。
> isDag.pl
#!/usr/bin/perl
use Graph;
use Getopt::Std;
getopt('i');
$file=$opt_i;
$g=Graph->new(%directed);
open(IN,"$file");
$n=0;
while($line=<IN>){
$line=~s/\n//;
($n1,$n2)=split("\t",$line);
$g->add_edge("$n1","$n2");
}
close(IN);
$isdag="ndag";
if($g->is_dag()){$isdag="dag";}
print "$isdag\n";
-----------------------
プログラムを実行してみる (実行にはGraphモジュールが必要です。CPANからgetしてください) 。
./isDag.pl -i example.txt
ndag
入力グラフは巡回 (あるノードを起点として同じノードへ戻ってくるパスがある) グラフなのでDAGではないので"ndag"と表示される（はず）。
>example.txt
nodeA<tab>nodeB
nodeB<tab>nodeC
nodeC<tab>nodeA
この２項関係をグラフで表すと"nodeA->nodeB->nodeC->nodeA"とcyclicなグラフとなります。このようなファイルを入力としてDAG (Directed Acyclic Graph : 非巡回有向グラフ)かどうかを判定するプログラムは以下のように簡単に書けます。
> isDag.pl
#!/usr/bin/perl
use Graph;
use Getopt::Std;
getopt('i');
$file=$opt_i;
$g=Graph->new(%directed);
open(IN,"$file");
$n=0;
while($line=<IN>){
$line=~s/\n//;
($n1,$n2)=split("\t",$line);
$g->add_edge("$n1","$n2");
}
close(IN);
$isdag="ndag";
if($g->is_dag()){$isdag="dag";}
print "$isdag\n";
-----------------------
プログラムを実行してみる (実行にはGraphモジュールが必要です。CPANからgetしてください) 。
./isDag.pl -i example.txt
ndag
入力グラフは巡回 (あるノードを起点として同じノードへ戻ってくるパスがある) グラフなのでDAGではないので"ndag"と表示される（はず）。
2008-05-18
偏相関係数
３変数の関係性を検討するとしましょう。例えば、以下のような因果関係があるとします。
我々はA,B,Cがなんなのかを知っていますし、観察や測定でそれらの値を測ることもできます。ただ、A,B,Cの間の因果関係については未知で、どうにかして測定データからそれを知りたいと考えます。こんなとき、変数間の関係性の強さの計量として相関係数がよく用いられます。Rで試してみます。
g<-function(n){
A<-rnorm(n,0,0.5)
B<-0.5*A+rnorm(n,0,0.1)
C<-0.5*B+rnorm(n,0,0.1)
m<-matrix(c(A,B,C),ncol=3,byrow=F)
m
}
d<-g(100)
100次元の人工データを生成しました、相関係数は以下のようになります。
cor(d)
A B C
A 1.0000000 0.9226645 0.9048547
B 0.9226645 1.0000000 0.8432231
C 0.9048547 0.8432231 1.0000000
BとCの相関係数は0.843です。独立性の検定を行えば、当然有意となります。図にあるように本当はBとCは独立なのにAを介して見せかけの相関が発生しています。これを偽相関といいます。これを回避する方法が偏相関係数です。偏相関係数は、他の変数値の影響を取り除いた２変数の直接相関を計算することができます。やってみましょう。
partial correlation
A B 0.697700 4.11E-17
A C 0.611920 3.93E-12
B C 0.050840 6.20E-01
自作のプログラムで計算しました。３列目が偏相関係数、４列めは独立性の検定におけるp値です。B,C間の相関が消失して正しい相関構造を捉えています。
偏相関係数の計算法はいくつかありますが、回帰を用いた方法を紹介します。
データカラムにA,B,Cが付加されているとします。まず、BをA, CをAで線形回帰します。
r1<-lm(d$B ~ d$A)
r2<-lm(d$C ~ d$A)
それぞれの回帰の残差の相関係数を計算します。
cor(r1$residuals,r2$residuals)
[1] 0.05083956
上記のBC間の偏相関係数と一致しました。つまり、２変数それぞれについて条件づける変数群で回帰し、残差間の相関係数が偏相関係数となります。
数千〜数万変数間の因果関係性を偏相関係数を基に高速に推定するソフトを開発しました。
我々はA,B,Cがなんなのかを知っていますし、観察や測定でそれらの値を測ることもできます。ただ、A,B,Cの間の因果関係については未知で、どうにかして測定データからそれを知りたいと考えます。こんなとき、変数間の関係性の強さの計量として相関係数がよく用いられます。Rで試してみます。
g<-function(n){
A<-rnorm(n,0,0.5)
B<-0.5*A+rnorm(n,0,0.1)
C<-0.5*B+rnorm(n,0,0.1)
m<-matrix(c(A,B,C),ncol=3,byrow=F)
m
}
d<-g(100)
100次元の人工データを生成しました、相関係数は以下のようになります。
cor(d)
A B C
A 1.0000000 0.9226645 0.9048547
B 0.9226645 1.0000000 0.8432231
C 0.9048547 0.8432231 1.0000000
BとCの相関係数は0.843です。独立性の検定を行えば、当然有意となります。図にあるように本当はBとCは独立なのにAを介して見せかけの相関が発生しています。これを偽相関といいます。これを回避する方法が偏相関係数です。偏相関係数は、他の変数値の影響を取り除いた２変数の直接相関を計算することができます。やってみましょう。
partial correlation
A B 0.697700 4.11E-17
A C 0.611920 3.93E-12
B C 0.050840 6.20E-01
自作のプログラムで計算しました。３列目が偏相関係数、４列めは独立性の検定におけるp値です。B,C間の相関が消失して正しい相関構造を捉えています。
偏相関係数の計算法はいくつかありますが、回帰を用いた方法を紹介します。
データカラムにA,B,Cが付加されているとします。まず、BをA, CをAで線形回帰します。
r1<-lm(d$B ~ d$A)
r2<-lm(d$C ~ d$A)
それぞれの回帰の残差の相関係数を計算します。
cor(r1$residuals,r2$residuals)
[1] 0.05083956
上記のBC間の偏相関係数と一致しました。つまり、２変数それぞれについて条件づける変数群で回帰し、残差間の相関係数が偏相関係数となります。
数千〜数万変数間の因果関係性を偏相関係数を基に高速に推定するソフトを開発しました。
2008-05-18
多重補正 (その３)
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
で参照できます。
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
で参照できます。
2008-05-18
しょぼ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とか)を使うときは注意。しかし、余計なことを...。
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とか)を使うときは注意。しかし、余計なことを...。
2008-05-18
Now available: MSigDB v2.5
2008-05-16
多重補正 (その２)
古典的なアプローチによる多重補正の計算は比較的簡単に計算できます。プログラミングも容易です。統計パッケージ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.
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.
2008-05-12
多重補正 (その１)
多群におけるペアワイズによる検定や独立要素の同一仮説の繰り返し検定は多重検定と呼ばれます。例えば、患者群A, B, Cそれぞれの平均血圧に差があるか？AとBでは？AとCでは？BとCでは？と繰り返すことは多重検定となります。また、1000個の遺伝子が解析対象で、群AとBにおける発現差が認められるか？といったとき遺伝子aは？遺伝子bは？... と繰り返します。これもまた多重検定。統計検定では有意水準αというものを設定します。よくあるのは0.01や0.05ですね。帰無仮説が偶然起こる確率です。つまり、この例の場合"遺伝子発現差がない"が起こる確率です。α=0.05をパスしたら5%有意といいます。逆にいうと5%は間違えるわけです。これが多重補正の概念のもと。繰り返し有意水準5%で検定を行うとリスクは5%ずつ蓄積していくという考え方です。よって、検定対象全体で有意性を語る場合、誤って有意！を避けるため、あらかじめ一回あたりの有意水準を検定回数で割っておこうという算段です。これが有名なBonferroniの補正です。理解しやすいように有意水準αは固定してp値の方を補正します 。
※nを検定回数、大文字Pを補正後のp値とします。
この手の補正後のpはFWER (Family Wise Error Rate) とも呼びます。しかし、こいつはあまりに厳しい補正ではありませんか。検定回数が1000回ならp=0.00005より小さくないと5%有意となりません。ゲノムワイドな解析でnが数万から数十万なんて場合は絶望的です。ただ、この方法はあまりに保守的です。補正が厳しすぎて本来は有意なものもみすみす逃してしまいます (偽陰性：false negative)。これらの補正法は各検定が独立という仮定をおいていて、独立でない場合、過度に補正をかけてしまいます。特に生物データ解析は各要素 (遺伝子等) がお互いに独立とは言い難いわけで、やりすぎちゃうわけです。また、"これ差があると思う"、"むしろ差が無いと困る"というのに補正で撃沈して論文にできない、補正しないとreviewerに怒られる、そんなの嫌だという大人の事情から (？) 、Bonferroni以降、多くの補正法が提案されています。
> Bonferroni
> Dunn-Sidak : Ury, H. K. (1976) A comparison of four procedures for multiple comparisons among means (pairwise contrasts) for arbitrary sample sizes, Technometrics, 18, 89-97.
Holm Holm, S. A. (1979) A simple sequentially rejective multiple test procedure, Scand J Stat, 6, 65-70.
> Simes : Simes, R. J. (1986) An improved Bonferroni procedure for multiple tests of significance, Biometrika, 73, 751-4.
> Hommel : Hommel, G. A. (1988) A stagewise rejective multiple test procedure based on a modified Bonferroni test, Biometrika, 75, 383-5.
> Hochberg : Hochberg, Y. (1988) A sharper Bonferroni procedure for multiple tests of significance, Biometrika, 75, 800-2.
> Benjamin & Hochberg : Benjamin, Y., Hochberg, Y. (1995) Controlling the False Discovery Rate: A practical and Powerful Approach to Multiple Testing, J R Stat Soc B, 57, 289-300.
これらの方法の詳細や背景はつづきで...
※nを検定回数、大文字Pを補正後のp値とします。
この手の補正後のpはFWER (Family Wise Error Rate) とも呼びます。しかし、こいつはあまりに厳しい補正ではありませんか。検定回数が1000回ならp=0.00005より小さくないと5%有意となりません。ゲノムワイドな解析でnが数万から数十万なんて場合は絶望的です。ただ、この方法はあまりに保守的です。補正が厳しすぎて本来は有意なものもみすみす逃してしまいます (偽陰性：false negative)。これらの補正法は各検定が独立という仮定をおいていて、独立でない場合、過度に補正をかけてしまいます。特に生物データ解析は各要素 (遺伝子等) がお互いに独立とは言い難いわけで、やりすぎちゃうわけです。また、"これ差があると思う"、"むしろ差が無いと困る"というのに補正で撃沈して論文にできない、補正しないとreviewerに怒られる、そんなの嫌だという大人の事情から (？) 、Bonferroni以降、多くの補正法が提案されています。
> Bonferroni
> Dunn-Sidak : Ury, H. K. (1976) A comparison of four procedures for multiple comparisons among means (pairwise contrasts) for arbitrary sample sizes, Technometrics, 18, 89-97.
Holm Holm, S. A. (1979) A simple sequentially rejective multiple test procedure, Scand J Stat, 6, 65-70.
> Simes : Simes, R. J. (1986) An improved Bonferroni procedure for multiple tests of significance, Biometrika, 73, 751-4.
> Hommel : Hommel, G. A. (1988) A stagewise rejective multiple test procedure based on a modified Bonferroni test, Biometrika, 75, 383-5.
> Hochberg : Hochberg, Y. (1988) A sharper Bonferroni procedure for multiple tests of significance, Biometrika, 75, 800-2.
> Benjamin & Hochberg : Benjamin, Y., Hochberg, Y. (1995) Controlling the False Discovery Rate: A practical and Powerful Approach to Multiple Testing, J R Stat Soc B, 57, 289-300.
これらの方法の詳細や背景はつづきで...
2008-05-11
Biomarker Discovery, Validation and Applications
今年の２月に国際シンポジウム"Biomarker Discovery, Validation and Applications"に参加してきました。バイオマーカー研究／開発に関係する欧米の大学・公的研究機関・企業等が集参加し、研究開発成果、バイオマーカー開発における問題点、規制側の取り組み等など、幅広い内容で構成されているシンポジウムでした。現在の主流はゲノムマーカーですが、シンポジウムではタンパク質マーカーに関する議論が主流でした。技術的な側面では特に革新的なトピックはなく、分析面では質量分析や抗体アレイが中心、統計・インフォマティクスは検定の多重補正等、退屈な内容でした。※ちなみに抗体アレイはまだまだ発展途上なイメージでした。ただ、製薬会社等の演題で、度々"Systems Biology"がキーワード的に使われていました。バイオマーカー開発においてもバイオマーカーの機序をシステムとして理解する必要があり、特に“がん”領域のバイオマーカーにおいては、シグナル伝達パスウェイを中心とした機序を踏まえたバイオマーカーの開発が今後の主流となる印象を受けました。創薬を意識したバイオマーカー開発においてテキストマイニングを積極的に用い、網羅解析によるデータを過去に報告された生物学的知見と統合し、システムとしての整合性を検証するような試みが行われており、予想以上に"Systems Biology"に対する意識が強い印象をうけました。今は、まだ地味なインフォマティクス技術が今後より重要になっていくと思います。
ところで、シンポジウムはアメリカのレイク・タホで行われました。有数のスキー・リゾート地だそうです。なので参加者の多くはスキーを楽しんでいたようです (スキーウェアで参加している人もいました笑) 。写真はレイク・タホの桟橋から撮ったものです。
帰りの飛行機が飛ばずサンフランシスコで足止めくらいました。いつも思いますがいい加減なシステムですね。日本のシステムの優秀さが再認識されます。SFでは山中で食べられなかったシーフードをScoma'sで食べました。※地元で有名な店のようです。いい雰囲気でした。写真はSFのフィッシャーマンズ・ワーフ。
ところで、シンポジウムはアメリカのレイク・タホで行われました。有数のスキー・リゾート地だそうです。なので参加者の多くはスキーを楽しんでいたようです (スキーウェアで参加している人もいました笑) 。写真はレイク・タホの桟橋から撮ったものです。
帰りの飛行機が飛ばずサンフランシスコで足止めくらいました。いつも思いますがいい加減なシステムですね。日本のシステムの優秀さが再認識されます。SFでは山中で食べられなかったシーフードをScoma'sで食べました。※地元で有名な店のようです。いい雰囲気でした。写真はSFのフィッシャーマンズ・ワーフ。
2008-05-10
KNIME
データ解析のパターンをテンプレート化するソフトウェアがあります。データの取り込みや書き出し、フォーマット変換やフィルタリング、統計解析や機械学習などデータ解析に必要な機能はノードとして具現化され、各ノード間をエッジで結ぶことにより解析フローを構築します。このようなフロー図をワークフローと呼びます。私の会社ではKNIME (Konstanz Information Miner)というワークフローソフトウェアのノード開発／販売を行っています。KNIMEは"ないむ"とよみます。ドイツのコンスタンツ大学が開発したフリーウェアでデータマイニング系の機能とケムインフォマティクスの解析機能が充実しています。バイオインフォマティクス関係の機能は無いのですが、JAVAを用いてノードを自作することができます。バイオマーカー開発にはRandom ForestやSVMなどの機械学習に関する機能が使えそうです。クロスバリデーションのようなバッチ機能もあり、RやWekaをハンドリングできるので解析機能のとしてはいい線いっていると思います。私は、普段スクリプトベースで仕事をしているので、ワークフローを直接使う動機付けがなかなかないのですが、確立した解析パターンやノウハウを共有するのに有効なツールだと思います。
2008-05-08
がんマーカーに求められる精度
前回、バイオマーカーの評価法として陽性予測値（PPV）について書きました。同じ感度・特異度でも対象疾患の事前確率（有病率）の違いで、臨床における実用性が大きく異なることを示しました。下表は主要ながんに対するスクリーニング目的のバイオマーカーに要求される特異度を計算したものです。
罹患率は国立がんセンター"がんの統計 ’07"を参照しました。罹患率とは同部位別がん罹患率です。この罹患率は10万人に対するものなので、1人あたりに換算したものを有病率としました。総じて低い有病率なので感度よりも特異度が要求されます (前回の検討参照) 。感度を0.9としたときの１．PPV=0.1即ち、10人の陽性判定のうち1ががんであるようなときの特異度、2．PPV=0.5、2人の陽性判定のうち 1ががんであるようなときの特異度を示しました。総じて、非常に高い特異度が要求されます。もっとも、有病率は専門病院においてはもっと高いことが推定されますが、あくまでも統計から大まかに計算を試みました。とはいえ、多様な背景を持つ対象群から高感度検出ねらうスクリーニングマーカーに要求される特異性はとてつもなく大きなものになると考えられます。次回は、実用化されている腫瘍マーカーについて検討してみたいと思います。
罹患率は国立がんセンター"がんの統計 ’07"を参照しました。罹患率とは同部位別がん罹患率です。この罹患率は10万人に対するものなので、1人あたりに換算したものを有病率としました。総じて低い有病率なので感度よりも特異度が要求されます (前回の検討参照) 。感度を0.9としたときの１．PPV=0.1即ち、10人の陽性判定のうち1ががんであるようなときの特異度、2．PPV=0.5、2人の陽性判定のうち 1ががんであるようなときの特異度を示しました。総じて、非常に高い特異度が要求されます。もっとも、有病率は専門病院においてはもっと高いことが推定されますが、あくまでも統計から大まかに計算を試みました。とはいえ、多様な背景を持つ対象群から高感度検出ねらうスクリーニングマーカーに要求される特異性はとてつもなく大きなものになると考えられます。次回は、実用化されている腫瘍マーカーについて検討してみたいと思います。
2008-05-05
解析ツール公開
このブログで取り上げた解析に関する自作ツールを公開しました。
まずは、クライアントサイド (Javascript) で計算可能なツールを実装しました。計算量の多いものはCGIで提供する予定です。逐次更新したいと思います。
まずは、クライアントサイド (Javascript) で計算可能なツールを実装しました。計算量の多いものはCGIで提供する予定です。逐次更新したいと思います。
2008-05-04
サイクリング
千住新橋から和光市まで荒川沿いのサイクリングロードを嫁と走ってきました。自宅から千住新橋までを加えると片道約23Kmです。
折り返した後、土手でおにぎりを食べて帰りました。心地よかったです。お尻が痛くなったのでパット付きパンツを買おうと思います。
我が愛車です。るいがの。
折り返した後、土手でおにぎりを食べて帰りました。心地よかったです。お尻が痛くなったのでパット付きパンツを買おうと思います。
我が愛車です。るいがの。
2008-05-02
バイオマーカーの切れ味
例えば、ある疾患を有する人と健常人を識別するようなスクリーニング目的のバイオマーカーの性能は以下のように評価します。
> TP (True Positive)：真陽性：病気の人を正しく陽性と判定
> TN (True Negative)：真陰性：健常人を正しく陰性と判定
> FP (False Positive)：偽陽性：健常人を誤って陽性と判定
> FN (False Negative)：偽陰性：病気の人を誤って陰性と判定
TP・TNがなるべく大きく、FP・FNがなるべく小さいマーカーが優秀といえます。この優秀さの指標として以下が使用されます。
> Sensitivity：感度：病気の人の中で陽性と判定される割合
> Specificity：特異度：健常人の中で陰性と判定される割合
どちらも０〜１を示し、ともに１のとき完全な判定が可能なマーカーといえます。
一般的に、バイオマーカーは臨床背景をコントールした上で限られた検体数を用いて開発されます。マーカーには高い感度と特異度が要求されますが、マーカーが適用される疾患の状況によっては、感度と特異度だけでマーカーの性能を議論することは危険です。以下の指標を評価する必要があります。
> PPV (Positive Predictive Value)：陽性予測値：検査が陽性のとき、実際に病気である割合
> NPV (Negative Predictive Value)：陰性予測値：検査が陰性のとき、実際に健常である割合
PPV は1に近いほど優れ、検査で陽性判定を受けた1/PPV人に1人が、実際に疾患を有することになります。このPPVには疾患に関する"事前確率"が多大な影響を与えます。スクリーニングの場合、この"事前確率"とは疾患の有病率となります。例えば感度が1・特異度が0.95の高性能なマーカーを開発したとします。感度が1なので病気の人を100%陽性判定し、健常人の95%を陰性判定します。このマーカーの適応疾患は有病率が1万人に1人（0.00001）の疾患だとします。仮に100万人（※何人でもいいのですが）にテストしたとすると結果は以下の表になります。
PPVが0.00020なので、1/0.0002=5001人に1が実際に疾患を有しています。よって、5001回テストしたとき5000回は、本当は病気でないのに陽性と判定されてしまいます。日々の検査では、ほとんどが偽陽性判定となるでしょう。
予測対象の疾患の有病率によって、感度と特異度の影響は異なります。下図は有病率0.01, 0.50, 0.99のときの感度（紫）・特異度（青）に対応する陽性予測値PPVの値です。
> 有病率：0.01
予測が難しいことがわかります。PPVは感度（縦軸）の影響は低いですが、特異度（横軸）に鋭敏です。元々、有病率が低いので、多くの健常人を正しく陰性判定するマーカーが要求されます。
> 有病率：0.5
有病率が50%の場合、PPVは感度・特異度ともに鋭敏です。
> 有病率：0.99
有病率が高い場合、PPVは特異度の影響は低いですが、感度に鋭敏です。既に判定対象は、ほぼ疾患を有しているので、その中から正しい陽性判定を下すマーカーが要求されます。
以上の結果をみる限り、レアな疾患をスクリーニングすることは (当然ですが) 容易でないことがわかります。がんの早期発見等を目論んで、血中マーカーの開発が行われていますが、血中タンパク質はダイナミックレンジが非常に大きく、微量のマーカータンパクの検出限界の問題が発生します (A) 。特殊な技術を用いてこれをクリアしたとして、次の問題は特異度です。いかに、対象がんを特異的に検出できるかです (B) 。これをクリアしたとして、最後に立ちはだかる問題が、上述した有病率です。レアながんの場合、いかなる高精度なマーカーを持ってしても、実際の検査で役に立たない可能性もあります (C) 。私が参加しているプロジェクトでは、術後（あるいは生検）組織に対して実施するマーカーを開発しています。腫瘍組織を対象としますので血中に比べると安定的に検査でき、感度も担保できる可能性が高いと考えます。※ (A) に対する優位性。 (B) 特異性については、血中マーカーに比べるとダイバーシティの次元が低く、はるかにハードルが低いと考えます。スクリーニングの特異性は、健常人や様々な疾患等の背景を持つ多様性に対して要求されるのに対して、がん組織マーカーの特異性は、予後の良／不良等のシンプルなものだからです (それでも十分難しいですが) 。 (C) に関しては、有病率という事前確率を1にできます。なぜなら、検査対象は腫瘍組織であり診断が確定済みだからです。要求されるのは高い感度です。これらの考察から、スクリーニング目的の血中マーカーは実現すれば、非常に強力ですが、実現は困難を極めます。一方、術後組織マーカーは、使用方法はフォーカスされていますが、実現可能性が高いと考えられます。いかがでしょうか？
※陽性・陰性予測値はこちらで計算できます。
> TP (True Positive)：真陽性：病気の人を正しく陽性と判定
> TN (True Negative)：真陰性：健常人を正しく陰性と判定
> FP (False Positive)：偽陽性：健常人を誤って陽性と判定
> FN (False Negative)：偽陰性：病気の人を誤って陰性と判定
TP・TNがなるべく大きく、FP・FNがなるべく小さいマーカーが優秀といえます。この優秀さの指標として以下が使用されます。
> Sensitivity：感度：病気の人の中で陽性と判定される割合
> Specificity：特異度：健常人の中で陰性と判定される割合
どちらも０〜１を示し、ともに１のとき完全な判定が可能なマーカーといえます。
一般的に、バイオマーカーは臨床背景をコントールした上で限られた検体数を用いて開発されます。マーカーには高い感度と特異度が要求されますが、マーカーが適用される疾患の状況によっては、感度と特異度だけでマーカーの性能を議論することは危険です。以下の指標を評価する必要があります。
> PPV (Positive Predictive Value)：陽性予測値：検査が陽性のとき、実際に病気である割合
> NPV (Negative Predictive Value)：陰性予測値：検査が陰性のとき、実際に健常である割合
PPV は1に近いほど優れ、検査で陽性判定を受けた1/PPV人に1人が、実際に疾患を有することになります。このPPVには疾患に関する"事前確率"が多大な影響を与えます。スクリーニングの場合、この"事前確率"とは疾患の有病率となります。例えば感度が1・特異度が0.95の高性能なマーカーを開発したとします。感度が1なので病気の人を100%陽性判定し、健常人の95%を陰性判定します。このマーカーの適応疾患は有病率が1万人に1人（0.00001）の疾患だとします。仮に100万人（※何人でもいいのですが）にテストしたとすると結果は以下の表になります。
PPVが0.00020なので、1/0.0002=5001人に1が実際に疾患を有しています。よって、5001回テストしたとき5000回は、本当は病気でないのに陽性と判定されてしまいます。日々の検査では、ほとんどが偽陽性判定となるでしょう。
予測対象の疾患の有病率によって、感度と特異度の影響は異なります。下図は有病率0.01, 0.50, 0.99のときの感度（紫）・特異度（青）に対応する陽性予測値PPVの値です。
> 有病率：0.01
予測が難しいことがわかります。PPVは感度（縦軸）の影響は低いですが、特異度（横軸）に鋭敏です。元々、有病率が低いので、多くの健常人を正しく陰性判定するマーカーが要求されます。
> 有病率：0.5
有病率が50%の場合、PPVは感度・特異度ともに鋭敏です。
> 有病率：0.99
有病率が高い場合、PPVは特異度の影響は低いですが、感度に鋭敏です。既に判定対象は、ほぼ疾患を有しているので、その中から正しい陽性判定を下すマーカーが要求されます。
以上の結果をみる限り、レアな疾患をスクリーニングすることは (当然ですが) 容易でないことがわかります。がんの早期発見等を目論んで、血中マーカーの開発が行われていますが、血中タンパク質はダイナミックレンジが非常に大きく、微量のマーカータンパクの検出限界の問題が発生します (A) 。特殊な技術を用いてこれをクリアしたとして、次の問題は特異度です。いかに、対象がんを特異的に検出できるかです (B) 。これをクリアしたとして、最後に立ちはだかる問題が、上述した有病率です。レアながんの場合、いかなる高精度なマーカーを持ってしても、実際の検査で役に立たない可能性もあります (C) 。私が参加しているプロジェクトでは、術後（あるいは生検）組織に対して実施するマーカーを開発しています。腫瘍組織を対象としますので血中に比べると安定的に検査でき、感度も担保できる可能性が高いと考えます。※ (A) に対する優位性。 (B) 特異性については、血中マーカーに比べるとダイバーシティの次元が低く、はるかにハードルが低いと考えます。スクリーニングの特異性は、健常人や様々な疾患等の背景を持つ多様性に対して要求されるのに対して、がん組織マーカーの特異性は、予後の良／不良等のシンプルなものだからです (それでも十分難しいですが) 。 (C) に関しては、有病率という事前確率を1にできます。なぜなら、検査対象は腫瘍組織であり診断が確定済みだからです。要求されるのは高い感度です。これらの考察から、スクリーニング目的の血中マーカーは実現すれば、非常に強力ですが、実現は困難を極めます。一方、術後組織マーカーは、使用方法はフォーカスされていますが、実現可能性が高いと考えられます。いかがでしょうか？
※陽性・陰性予測値はこちらで計算できます。
2008-05-01
マイクロアレイデータ統計解析プロトコール
マイクロアレイデータ統計解析プロトコール 〜 Excelを中心としたデータの標準化から有意差解析，クラスタリング，ネットワーク解析法のすべて. 羊土社
が出版されました。マイクロアレイデータの解析に関して、前処理、発現差解析、クラスター解析、ネットワーク解析が網羅されています。特にネットワーク解析に関する実践的な解説書としては、他に無い内容となっています。私も部分的ですが"応用解析編：ネットワークと実験値の適合度評価"と題して、文献ベースの相互作用ネットワークと対応する測定データから整合性を統計的に評価するトピックについて書かせていただきました。よろしければ一読を。
