2008-06-05
偏相関係数（その２）
以前、偏粗関係数をもとに条件付き独立性を評価することにより、相関係数では検出されてしまう偽相関を取り除くことができることをしめしました。偏相関係数をみることで直接的な強い相関構造を得る事ができます。しかし、偏相関係数による直接相関の検出は必要条件ではあるが十分条件ではありません。
図のような因果関係をもとにデータを発生させてみます。
g<-function(n){
A<-rnorm(n,0,0.5)
B<-rnorm(n,0,0.5)
C<-0.5*A+0.5*B+rnorm(n,0,0.1)
m<-matrix(c(A,B,C),ncol=3,byrow=F)
m
}
d<-g(100)
cor(d)
[,1] [,2] [,3]
[1,] 1.0000000 0.0419679 0.6797696
[2,] 0.0419679 1.0000000 0.7144819
[3,] 0.6797696 0.7144819 1.0000000
相関係数行列により変数１と２は独立ということがわかります。では、変数３を条件付けた変数１と２の偏相関係数を計算してみます。
r1<-lm(d[,1] ~ d[,3])
r2<-lm(d[,2] ~ d[,3])
cor(r1$residuals,r2$residuals)
[1] -0.8646989
相関係数0.0419679であった変数１と２の間の偏相関係数は -0.8646989と強い負の関係性が生じてしまっています。これを”選択の偏り”と呼び、図のような合流構造で合流点を固定すると起こる現象です。例え、多次元正規分布しているデータでも偏相関係数だけみていれば直接相関構造が担保されるわけではないので注意です。
図のような因果関係をもとにデータを発生させてみます。
g<-function(n){
A<-rnorm(n,0,0.5)
B<-rnorm(n,0,0.5)
C<-0.5*A+0.5*B+rnorm(n,0,0.1)
m<-matrix(c(A,B,C),ncol=3,byrow=F)
m
}
d<-g(100)
cor(d)
[,1] [,2] [,3]
[1,] 1.0000000 0.0419679 0.6797696
[2,] 0.0419679 1.0000000 0.7144819
[3,] 0.6797696 0.7144819 1.0000000
相関係数行列により変数１と２は独立ということがわかります。では、変数３を条件付けた変数１と２の偏相関係数を計算してみます。
r1<-lm(d[,1] ~ d[,3])
r2<-lm(d[,2] ~ d[,3])
cor(r1$residuals,r2$residuals)
[1] -0.8646989
相関係数0.0419679であった変数１と２の間の偏相関係数は -0.8646989と強い負の関係性が生じてしまっています。これを”選択の偏り”と呼び、図のような合流構造で合流点を固定すると起こる現象です。例え、多次元正規分布しているデータでも偏相関係数だけみていれば直接相関構造が担保されるわけではないので注意です。
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-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-04-26
d-sep test
前回、graphical modelについて書きましたが、GMにより因果関係をグラフとして表すことができます。AがBの原因だと考えられれば A→B と表すのです。この二項関係を積み重ねていくとグラフになります。
さて、今、いくつかの観測（測定）可能な対象があって、それらの因果関係を知りたいとします。それらの対象（変数）の背景や事前知識から変数間の因果関係をグラフで表現します。このグラフが仮説となります。仮説を検証するために各変数のデータをとります。例えば、３遺伝子A, B, C間の転写制御関係を考えます。各遺伝子に対する事前知識よりA→B, A→Cといった仮説をたてることができます。３遺伝子の発現を測定しデータを得ます。この転写制御の仮説は正しいか統計的に検証しようという具合です。
このような状態でグラフが表す因果関係が正しいかどうか統計的に評価する方法は種々あります。その１つに"d-sep test"があります。d-sep testはグラフ構造上における条件付き独立性を検定する手法です。アルゴリズムは以下の通りです。
1. エッジがないノードペア i, j について
2. iとjのip, jpに対する条件付き独立性を検定しp値を算出する (ipはノードiのすべての親ノード)
3. 1.を満たすすべてのノードペアについて2.によりp値を算出する
4. 算出されたp値の個数をNとする
5. すべてのp値の対数をとり、その総和をとり-2倍する (C=-2Σlog(p))
6. Cを自由度２Nでカイ二乗検定しp値を得る
※3.のp値の算出法は変数が連続変数か離散変数かで異なる。連続変数の場合、i,j,ip,jpで偏相関係数行列を算出し、i-j要素をZ変換し正規分布で検定する。離散変数の場合、ip,jpで層別したi-jの分割表をカイ二乗検定する。
※5.はメタアナリシスのテクニックでFisher Cと呼ばれる複数のp値を統合する方法。
d-sep testは比較的簡単な検定により高速にグラフの整合性 (graph consistency) をチェックできる方法です。この手法はB Shipleyにより提案され、生物統計データに活用されています。
詳細は以下の書籍を参照ください。
Cause and Correlation in Biology: A User's Guide to Path Analysis, Structural Equations and Causal Inference
本書籍は生物データ解析における関連・因果解析の方法論を丁寧に説明してある良書です。具体例が多く、わかりやすく、即実践できます。是非一読を。
さて、今、いくつかの観測（測定）可能な対象があって、それらの因果関係を知りたいとします。それらの対象（変数）の背景や事前知識から変数間の因果関係をグラフで表現します。このグラフが仮説となります。仮説を検証するために各変数のデータをとります。例えば、３遺伝子A, B, C間の転写制御関係を考えます。各遺伝子に対する事前知識よりA→B, A→Cといった仮説をたてることができます。３遺伝子の発現を測定しデータを得ます。この転写制御の仮説は正しいか統計的に検証しようという具合です。
このような状態でグラフが表す因果関係が正しいかどうか統計的に評価する方法は種々あります。その１つに"d-sep test"があります。d-sep testはグラフ構造上における条件付き独立性を検定する手法です。アルゴリズムは以下の通りです。
1. エッジがないノードペア i, j について
2. iとjのip, jpに対する条件付き独立性を検定しp値を算出する (ipはノードiのすべての親ノード)
3. 1.を満たすすべてのノードペアについて2.によりp値を算出する
4. 算出されたp値の個数をNとする
5. すべてのp値の対数をとり、その総和をとり-2倍する (C=-2Σlog(p))
6. Cを自由度２Nでカイ二乗検定しp値を得る
※3.のp値の算出法は変数が連続変数か離散変数かで異なる。連続変数の場合、i,j,ip,jpで偏相関係数行列を算出し、i-j要素をZ変換し正規分布で検定する。離散変数の場合、ip,jpで層別したi-jの分割表をカイ二乗検定する。
※5.はメタアナリシスのテクニックでFisher Cと呼ばれる複数のp値を統合する方法。
d-sep testは比較的簡単な検定により高速にグラフの整合性 (graph consistency) をチェックできる方法です。この手法はB Shipleyにより提案され、生物統計データに活用されています。
詳細は以下の書籍を参照ください。
Cause and Correlation in Biology: A User's Guide to Path Analysis, Structural Equations and Causal Inference
本書籍は生物データ解析における関連・因果解析の方法論を丁寧に説明してある良書です。具体例が多く、わかりやすく、即実践できます。是非一読を。
