[連載] フリーソフトによるデータ解析・マイニング 第17

Rと判別分析

 


1.判別分析とは

私たち人間は毎日五感を通じて入力される膨大なデータを処理している。その中で最も多いのは、識別(discrimination)、分類(classification)、認識(recognition)に関する処理である。例えば、新聞や本などを読むときには、視覚を通じて入力されたデータと学習したデータとの照合を行い、その文字の読み方、文字・単語の意味などを識別・認識する。識別に関する能力は人間のみならず、他の動物も持っている。

このような識別・認識に関することを機械的に実現する研究分野がパターン認識(pattern recognition)である。パターン認識の典型的な例としては、郵便番号による手紙の自動分類や指紋・顔照合によるセキュリティ管理などがあげられる。

パターン認識は、コンピュータに事前に入力・記憶させたデータと識別すべきデータとの一致度を何らかのモデルによって計算する。その中、最も古典的な手法は、判別分析(discriminant analysis)である。

判別分析は、個体(あるいは対象が)がどのグループに属するかが明確である学習データを用いて判別モデルを構築し、そのモデルを用いて所属不明の個体(テスト用データ)がどのグループに帰属するかを判別する方法である。

古典的判別分析には、距離(あるいは類似度)による判別と判別関数による判別分析などの方法がある。

判別分析では、所属不明の個体が2つのグループのいずれに属するかを判別する問題を2群判別分析、3つ以上のグループのいずれかに属するかに関する判別問題を多群判別分析とよぶ。

判別問題における学習データは、次の表のようにどの個体がどのグループに属するかに関する情報が必要である。グループに関する情報は、回帰分析の被説明変数に対応するが、回帰分析の場合は被説明変数が量的であるのに対して、判別分析のグループに関する変数は質的変数である。

 

1 グループに関する学習用データ

個体

          

グループ情報

 

 

 表1の中のはグループに属する番目の個体の番目の変数である。グループに属する個体は個、グループに属する個体は個である。

 判別分析の問題では、グループを識別するためのモデル(規則、あるいは関数など)を求めるためのデータセット(学習データ)とそのモデルを用いて判別、あるいは評価を行うデータセット(テストデータ)が必要である。

 

2.距離による判別分析

ここでは、説明の便利のため2群判別分析を例とする。距離による判別分析では、まず2つのグループの中心を求める。中心として最も多く用いられているのはグループの平均ベクトルである。

グループのそれぞれの平均ベクトルをとする。

ここのは次のように定義されている。

(ただし=1〜)

 

次に所属不明の個体からグループの中心までの距離を求め、所属不明の個体は距離の小さい方のグループに属すると判断する。

今、所属不明の個体があるとする。

個体からグループの中心までの距離を、グループの中心までの距離をとする。

距離による判別分析では、ならば個体はグループに、ならば個体はグループに属すると判断する。

判別すべき個体が2つのグループの中心から等距離にある特殊な場合は、判別不能である。

 このような距離を用いた判別分析法は、グループ数が3以上の場合にも簡単に拡張することができる。またデータに関しては、どのような確率分布に従っているかのような条件を必要としないのが長所である。

距離の測度としては、Rにはユークリッド距離、市街距離、マハラノビス距離などの関数が用意されている。

マハラノビス(mahalanobis)距離は,多変量データ解析の書物に必ず登場する距離で,マハラノビス汎距離とも呼ぶ。マハラノビス距離は次ぎのように定義されている。

式の中の*はデータセット、はグループの中心、はそのグループの分散共分散行列、は分散共分散行列の逆行列である。

 Rにはマハラノビス距離を求める関数mahalanobisがある。その書き式を次に示す。

 

mahalanobis(,,)

 

実例を用いて説明するため、Rの中に用意されているirisデータを用いることにする。データirisに関しては本誌のNo.117(p.87)に説明されている。help(iris)コマンドを実行すると英文による説明を読むことができる。

ここでは問題を簡単にするため、irisデータの第1行から第50行までのsetosaという品種、第101行から第150行までのvirginicaという品種のデータを用いる。

  データirisの第1列から第4列まではそれぞれアヤメのがく片の長さと幅、花弁の長さと幅に関する計測データで、第5(Species)は品種(グループ)を表す質的データである。

 

>data(iris)#バージョン2.0以後は必要ではない

>iris[c(1,51,101),]

    Sepal.Length Sepal.Width Petal.Length Petal.Width    Species

1            5.1         3.5          1.4         0.2     setosa

51           7.0         3.2          4.7         1.4 versicolor

101          6.3         3.3          6.0         2.5  virginica

 

まずirisの中のsetosavirginicaの品種別のデータセットsetovirgiを作成する。品種ごとにデータセットを作成するので、グループの属性を示すSpecies列は必要がない。

 

> seto<-subset(iris[1:50,],select=-Species)

> virgi<-subset(iris[101:150,],select=-Species)

 

setovirgiはそれぞれ50行のデータである。ここでは、各データセットの145行を学習データとし、残りの5(4650)をテストデータとする。

次に学習用データの平均ベクトルと分散共分散行列を求める。

 

> seto.m<-apply(seto[1:45,],2,mean)

> virgi.m<-apply(virgi[1:45,],2,mean)

>seto.v<-var(seto[1:45,])

>virgi.v<-var(virgi[1:45,],)

                  

次のコマンドを実行するとテストデータseto[46:50]と両グループ(seto[1:45,]virgi[1:45,])の平均とのマハラノビス距離が求められる。

 

>D1<-mahalanobis(seto[46:50,],seto.m,seto.v)

>D2<-mahalanobis(seto[46:50,],virgi.m,virgi.v)

 

 求めた距離からseto[46:50,]setovirgi

の両グループの中のどのグループに属するかを確認するため、両グループまでの距離を次に示す。

 

> cbind(D1,D2)

        D1       D2

46 2.1752192 137.9376

47 2.8163645 173.8815

48 1.4346178 142.1425

49 1.2398930 182.5972

50 0.4700029 160.2070

 

 D1seto[46:50,]setoの中心、D2seto[46:50,]virgiの中心との距離である。各行のD1列の値がD2列の値より小さいのでseto[46:50,]setoのグループに属すると判断する。

同様な方法でvirgi[46:50,]と両グループの中心までのマハラノビス距離を求め、どのグループに属するかを判別することができる。

 説明の便利のため、ここでは学習用のデータとテスト用のデータを恣意的に前の45行と後の5行に分けているが、現実問題では、もっと説得力のある方法を取るべきである。

 マハラノビス距離の定義からわかるように、マハラノビス距離を求めるためには、学習データの分散共分散行列の逆行列を求める必要がある。しかし、すべてのデータの分散共分散行列の逆行列が求められるとは限らない。

 

3.判別関数による判別分析

判別関数による判別分析は、線形関数と非線形関数による判別分析に分けられる。

 

(1) 線形判別分析

線形判別分析(Linear Discriminant Analysis)は、グループ分けの境界が直線、あるいは超直面であり、次のような線形関数を用いてグループの所属の判別を行う方法である。

 

判別関数=

 

判別関数は、グループ間の分散とグループ内の分散の比を最大化することで求めることができる。

RのパッケージMASSの中にこのアプローチによる判別分析の関数オブジェクトldaがある。

 

>library(MASS)

 

関数ldaの詳細については、help(lda)で読むことができる。次にldaの最も簡潔な書式を示す。

 

   lda(formula, data)

 

 formulaには、回帰分析の場合と同様、「グループの識別変数~ 変数」のように記述する。

 ここでもRの中のirisデータを用いて、ldaの使用法について説明する。

 

@ データセットの作成

 前項では説明の便利のため、irisの2種類のデータを用いたがここでは3種類のirisを用いて3群判別分析を行う。

 まずデータirisから、学習データとテストデータを作成する。その作成方法はいろいろ考えられるが、ここでは奇数行と偶数行に二分することにする。

 

 

>even.n<-2*(1:75)-1

>train.data<-iris[even.n,]

>test.data<-iris[-even.n,]

 

上記のコマンドで150行の中から行番号が奇数である75行を学習データ(train.data)、残りの偶数行のデータをテストデータ(test.data)とするデータセットが作成される。

この段階で判別分析を行っても良いが、データirisの中のグループを示すラベルが若干長いので、散布図を作成すると見苦しくなる。そこで、setosaSversicolorCvirginicaVで示すことにする。

 

>Iris.lab<-factor(c(rep("S",25),rep("C",25),rep("V",25)))

>train.data[,5]<-Iris.lab

> train.data[c(1,26,51),]

   Sepal.Length Sepal.Width Petal.Length Petal.Width Species

1            5.1         3.5          1.4         0.2       S

51           7.0         3.2          4.7         1.4       C

101          6.3         3.3          6.0         2.5       V

 

A 判別関数を求める

学習データセットを次ぎのように用いると判別分析に必要となる統計量が求められる。ただし、次の書き式は学習データセットの中の変数をすべて用いた場合の書き方である。学習データセットの中の変数の1部分を用いる際には、のように用いる変数を記号「+」でつなぐ。

 

>(Z<- lda(Species~ .,data=train.data))

Call:

lda(Species ~ ., data = train.data)

 

Prior probabilities of groups:

        C         S         V

0.3333333 0.3333333 0.3333333

 

Group means:

  Sepal.Length Sepal.Width Petal.Length Petal.Width

C        5.992       2.776        4.308       1.352

S        5.024       3.480        1.456       0.228

V        6.504       2.936        5.564       2.076

 

Coefficients of linear discriminants:

                    LD1        LD2

Sepal.Length -0.5917846 -0.1971830

Sepal.Width  -1.8415262  2.2903417

Petal.Length  1.6530521 -0.7406709

Petal.Width   3.5634683  2.6365924

 

Proportion of trace:

   LD1    LD2

0.9913 0.0087

 

返された線形判別係数(Coefficients of linear discriminants)を用いて判別関数を構築する。判別関数は、用いる変数と判別係数との線形結合である。

この問題では2組の判別係数(LD1、LD2)が返されている。この判別係数を用いた第1判別関数を次に示す。

 

=

 

判別関数式の中のはそれぞれirisデータの変数Sepal.Length、Sepal.Width、 Petal.Length、Petal.Widthを示す。

定数項はグループの平均と判別係数を次ぎのように用いて求めることができる。第1列の値は第1判別関数の定数項で、第2列は第2判別関数の定数項である。

 

> apply(Z$means%*%Z$scaling,2,mean)

     LD1      LD2

1.486146 6.282412

 

上記の判別関数を用いた判別得点はコマンドpredict(Z)$xで返される。

判別関数で得られた値を判別得点(discriminant score)とよぶ。各判別関数が全体のグループ間分散をどのぐらい説明できるかはグループ間分散の比率から読み取れる。

 

B 学習データにおける判別結果

判別関数を用いて、学習データについて判別を行った結果は関数predictを用いて返す。関数predictは次の値を返す。

 

predict()$class

predict()$posterior

predict()$x

 

$classは各個体が判別されたグループのラベルで、$posteriorは各個体がどのグループに判別されているかに関する事後確率(0〜1)、$xは各個体の判別関数得点である。

学習データにおける判別結果は次のような表で確認することができる。

 

> table(train.data[,5],predict(Z)$class)

    C  S  V

  C 24  0  1

  S  0 25  0

  V  1  0 24

 

Cグループの1つがVグループに、Vグループの1つがCグループに誤判別され、誤判別率は2/75≒0.0267(2.67%)である。どの個体が誤判別されているかは、次のようなコマンドで追跡することができる。

 

>data.frame(train.data[,5],predict(Z)$class)

 

判別関数得点をグラフに示すこともできる。次のコマンドで第1判別関数得点のグループごとのヒストグラム(分布)が作成される。

 

>plot(Z,dimen=1)

 

作成されたヒストグラムから分かるように、S(setosa)はC(versicolor)、V(virginica)と重ならないが、C(versicolor)の右辺とV(virginica)

の左辺は若干重なる。重なる領域が多いほどお互いに間違って判別される確率(誤判別率)が高い。

 

1 学習データの第1判別関数得点の分布

 

  引数dimen=2を用いると、図2のような、横軸を第1判別関数、縦軸を第2判別関数とした散布図が作成される。散布図からグループ間の分類状況がマクロ的に把握できる。判別関数が3つ以上の場合は、dimenの値を3以上にすると対散布図が作成される。

 

>plot(Z,dimen=2)

2 学習データの判別関数得点の散布図

 

C 判別関数による判別

学習データで求めた判別関数を用いてテストデータについての判別分析は次ぎのように関数predictを用いる。

 

>Y<-predict(Z,test.data)

 

どれぐらい正しく判別されているかは、テストに用いたデータのグループ情報と判別結果のグループ情報のクロス表を作成することで確認できる。

 

>table(test.data[,5],Y$class)       

    C  S  V

  C 24  0  1

  S  0 25  0

  V  2  0 23

 

 上記の結果から分かるように、テストデータの中のS(setosa)はすべて正しく判別され、C(versicolor)は1つがV(virginica)と誤判別され、V2つがCと誤判別されている。誤判別されているのは3個で全体の中で占める割合(誤判別率)3/750.04(4)、正判別率は1-0.04=0.96 (96%)である。

テストデータの判別得点の散布図を用いて、判別状況を視覚的に考察することもできる。

 

>plot(Y$x,type="n")

>text(Y$x,labels=as.character(Iris.lab))

 

3 テストデータの判別関数値の散布図

 

 このアプローチによる線形判別分析は、各グループの母分散が等しいであるとの仮定に基づいている。各グループの母分散が異なる場合は、マハラノビス距離による判別のような母分散に制約がない方法を用いるべきである。

 

(2) 非線形判別分析

 古典的な非線形判別分析では、2次式を含む非線形関数を用いる方法がある。RのパッケージMASSの中には2次式による判別分析の関数qdaが用意されている。 関数qdaの使用方法は、ldaほとんど同じである。

 一般化線形モデルで紹介した、目的変数が質的変数であるロジスティック回帰も非線形判別分析の一つの方法である。

 判別関数のアプローチによる判別分析は、属性を示す変数の数が多くなると判別関数の構築が難しくなる。

 

. 交差確認

 判別分析やパターン認識の分野では、データセットから学習用のデータとテスト用のデータに分けてモデルの構築と評価をするのに交差確認(cross validation:交差検証、交差妥当化とも訳されている)という方法が広く知られている。

 交差確認では、データセット全体をn部分(サブデータセット)に均等に分割し、そのうちの1つをテスト用のデータとして残し、それ以外のn-1つを学習用のデータとして用いる。 

データセットをn部分に分割したときをn重交差確認(n-fold cross validation)法と呼ぶ。n重交差確認法では、一つのデータセットに対し、n回のモデルの構築とテスト(確認、検証)を行い、そのn回のテスト結果を全体の評価に用いる。

関数ldaqdaでは、データセットから1つの個体を除いて学習を行い、学習データに用いていない1つの個体で判別モデルの評価を行う作業を、すべての個体に対して繰り返す交差確認(leave-one-out cross-validation)の引数CVが用意されている。これは、n重交差確認のnが個体の数に等しい特殊なケースである。

デフォルトでCV=FALSEになっている。引数CV=Tにすると1つを除いた交差確認による結果が返される。次にirisデータを用いた1つを除いた交差確認の結果を示す。

 

> iris.CV<-lda(Species~.,data=iris,CV=T)

>(lda.tab<- table(iris[,5],iris.CV$class))

      setosa versicolor virginica

setosa     50      0          0      

versicolor  0     48          2      

virginica   0      1         49

 

 

判別率と誤判別率は次のコマンドで求めることができる。

 

> sum(lda.tab[row(lda.tab)==col(lda.tab)])/sum(lda.tab)

[1] 0.98

> sum(lda.tab[row(lda.tab)!=col(lda.tab)])/sum(lda.tab)

[1] 0.02

 

 個体数が十分大きくない古典的な判別分析では、1つの個体を除いた交差確認法が多く用いられているが、ニューラルネットワーク、決定木、サポートベクトルマシンのようなパターン認識の方法では、一般のn重交差確認法が多く用いられている。

nをいくつにするべきであるかは、用いるデータのサイズに依存するため明確な基準がなく、n=2、3、4、5、10が多く用いられているのが現状である。

残念ながら、判別関数lda、qdaにはleave-one-out cross-validation以外の交差確認の機能が用意されていない。