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.
