まめトーーク！

家の近くのショッピングモールにイアン・ソープが来ました。ヤクルトのスポーツ飲料のプロモーションで。事前に告知していなかったようで、人はそれほど多く集まっていなかった。トークショーとクイズ大会と握手会でした。握手後、サイン色紙を貰いました。すでに引退してますが、体つきが水泳に最適化されている印象。胸板がとてもあつく、なで肩、特に印象的だったのがひざです。ひざ間接が柔らかいので、逆（前）方向に大きく曲がっていました。２６歳だそうです。好青年といった感じでした。

以前、偏粗関係数をもとに条件付き独立性を評価することにより、相関係数では検出されてしまう偽相関を取り除くことができることをしめしました。偏相関係数をみることで直接的な強い相関構造を得る事ができます。しかし、偏相関係数による直接相関の検出は必要条件ではあるが十分条件ではありません。



図のような因果関係をもとにデータを発生させてみます。

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と強い負の関係性が生じてしまっています。これを”選択の偏り”と呼び、図のような合流構造で合流点を固定すると起こる現象です。例え、多次元正規分布しているデータでも偏相関係数だけみていれば直接相関構造が担保されるわけではないので注意です。

連続値をとる変数間の (条件付き) 独立性の指標は (偏) 相関係数です (ただし多次元正規分布の仮定のもと)。偏相関係数の計算／検定法は割とよく知られているので、ここでは離散値 (カテゴリカル) をとる変数間の (条件付き) 独立性の検定法を紹介します。
帰無仮説が真のとき、以下のG2統計量はカイ2乗分布に従います。




ここで、N_ijkはセルijkの観測頻度、m_ijkはセルijkの期待頻度です。kはSのカテゴリーが取り得るすべてのパターンをなめます。自由度dfは(X_iのカテゴリー数-1)*(X_jのカテゴリー数-1)とSのすべての要素のカテゴリー数の積です。このG2統計量を用いて周辺・条件付独立性を検定するには、自由度dfのカイ二乗分布におけるパーセンタイルを考えればよいことになります。



のときX_iとX_jは独立と判定します。



のときX_iとX_jはSに対して条件付き独立と判定します。

十〜数百、それ以上のデータがある状態で分類器を構築する場合、私は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)

グラフ構造に対するアルゴリズムを実装するには、グラフを表現するデータ構造を選択する必要があります。ノード間の関係は行列によりあらわすことができます。これは隣接行列と呼ばれ、例えば、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"と表示される（はず）。