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

Rと対応分析

1.対応分析とは

  対応分析 (correspondence analysis) は、フランスのベンゼクリ (Benzecri) によって1960年代に提唱され、1970年代から普及し始めたカテゴリカルデータの解析方法で、コレスポンデンス分析とも呼ばれている。
  類似の方法としては、1940年代に林知己夫氏によって提案された数量化V類、1980年代に西里静彦氏によって提案された双対尺度法 (dual scaling) などがある。それぞれの方法が提案された背景は異なるが、基本的なアプローチおよびアルゴリズムの中核は同じである。
 データ形式によっては、それぞれの手法の解析結果は変換によって一致させることも可能である。一時的には、数量化V類と対応分析は異なるデータ分析方法と見なされたが、既に数理的には同等であることが証明されている。
  数量化V類および対応分析の基本的考え方は、分割表において、行の項目と列の項目の相関が最大になるように、行と列の双方を並び替えることである。問題解決のアプローチは、主成分分析、因子分析とほぼ同じである。そこで、対応分析を制約つき主成分分析 [1]、正準相関分析 [2]、質的因子分析として見なす研究者もいる。
  対応分析は、データの構造を再現する面では、古典的主成分分析より効果が劣るが、特徴別に分類する面では、古典的主成分析より良い結果を示すケースが多い。
  対応分析の大まかなアルゴリズムを説明するため、表1のような度数に関する分割表があるとする

表1 データ行列F r×c
x1 x2 xj xc 合計
個体1 f11 f12 f1j f1c f1●
個体2 f21 f22 f2j f2c f2●
個体i fi1 fi2 fij fic fi●
個体r fr1 fr2 frj frc fr●
合計 f●1 f●2 f●j f●c n

 データ行列F の各要素を総度数n で割ったデータ(P )を表2に示す。

Pij=fij/n

表2 データ行列 Pr×c
x1 x2 xj xc 合計
個体1 p11 p12 p1j p1c p1●
個体2 p21 p22 p2j p2c p2●
個体i pi1 pi2 pij pic pi●
個体r pr1 pr2 prj prc pr●
合計 p●1 p●2 p●j p●c n

  このような分割表の独立性の問題では、通常カイ2乗統計量を用いる。カイ2乗統計量の各セルの値は次の式で定義されている。

カイ2乗統計量の各セルの値の式

カイ2乗統計量の各セルの値の式

  対応分析では、データ Fr×c あるいは P r×c を次の式で変換したデータの固有値問題に帰する。

Zijの式

表3 データ行列 Zr×c
x1 x2 xj xc 合計
個体1 z11 z12 z1j z1c z1●
個体2 z21 z22 z2j z2c z2●
個体i zi1 zi2 zij zic zi●
個体r zr1 zr2 zrj zrc zr●
合計 z●1 z●2 z●j z●c n

  対応分析では、分割表の列の効果はQ = Z t Z、行の効果はQ t = Z Z t の固有ベクトル(U,V ) をそれぞれ で基準化した値 ( D c-1/2UD r-1/2V ) を用いる。ここの Dcp●j を要素とした対角行列、D rpi● を要素とした対角行列である。
  Q = Z tZ の固有値が λ 1 ≧ λ 2 ≧ … ≧ λ kk = rank(Q ) ≦ min ( r , c ) である場合、次の式が成り立つ。

x2乗の式

2.パッケージMASSの対応分析

  2元分割表のデータを対象とする対応分析をシンプルコレスポンデンス (Simple Correspondence) 分析、多項目の反応パターン表のデータを対象とした対応分析をマルチプルコレスポンデンス (Multiple Correspondence) 分析と呼ぶが、通常前者を略して対応分析、後者を多重対応分析と呼ぶ。
  Rには、幾つかのパッケージに対応分析に関連する関数が含まれているが、まずパッケージMASSの中の関数 corresp、mca を説明することにする。

2.1 対応分析

  パッケージMASSの中の対応分析の関数 corresp の書き式を次に示す。

corresp(x,nf,…)

  引数xはデータマトリックスあるいはデータフレームである。nf は、求める軸の数(主成分、因子とも呼ぶ)を指定する引数である。デフォルトは1になっているのでnfを省略した場合は、1軸の結果のみが返される。
  関数 corresp を用いるためにはパッケージMASSを読み込む (library(MASS)) 手続きが必要である。
  データを用いて関数 corresp の使用方法を説明することにする。データの入力の手間を省くため、パッケージMASSの中のデータセット caith を用いることにする。
  データ caith はイギリスに住んでいる人々の目の色 (blue, light, medium, dark) と髪の色に (fair, red, medium, dark, black) 関して5387人を対象として行った調査結果である。
  データセット caith は4行5列の2元分割表である。行が目の色、列が髪の色になっている。

> caith

fair red medium dark black
blue
326 38 241 110 3
light
688 116 584 188 4
medium
343 84 909 412 26
dark
98 48 403 681 85

  関数 corresp を用いたデータ caith の対応分析のコマンドおよびその出力結果を次に示す。

>(caith.ca<-corresp(caith,nf=4))

First canonical correlation(s): 4.463684e-01 1.734554e-01 2.931691e-02 1.134031e-16

Row scores:

[,1] [,2] [,3]   [,4]
blue
-0.89679252 0.9536227 2.1884132 1
light
-0.98731818 0.5100045 -1.0837859 1
medium
0.07530627 -1.4124778 0.1894089 1
dark
1.57434710 0.7720361 -0.1482208 1

Column scores:

[,1] [,2] [,3] [,4]
fair
-1.21871379 1.0022432 0.4271282 -0.8692696
red
-0.52257500 0.2783364 -4.0268545 -1.3400421
medium
-0.09414671 -1.2009094 0.1103959 -0.8453208
dark
1.31888486 0.5992920 0.3450676 -1.2251588
black
2.45176017 1.6513565 -1.5736976 1.1609621

  ここでは軸数の引数nfを4にしている。関数 corresp では、累積寄与率を返さないので累積寄与率を計算するためには、nf=min (行数、列数)にしたほうがよい。もしデータ行列のランクが min (行数、列数)より小さいときには、nf を下げればよい。
  返される結果は、正準相関 (canonical correlation)、行の得点 (Row scores)、列の得点 (Column scores) の順になっている。返される正準相関は大きい順に並べられ、その2乗が固有値に等しい。対応分析では、計算された軸の行・列に対応する値をそれぞれ行の得点、列の得点と呼ぶ。
  対応分析でも主成分分析や因子分析と同じく各軸の寄与率や累積寄与率について考察を行うが、関数 corresp は寄与率を返さないので、正準相関を用いて算出する必要がある。
  正準相関は $cor に記録されている。データ caith の固有値および寄与率は次のコマンドで求めることができる。

> 固有値<-caith.ca$cor^2 
> round(固有値,3)
[1] 0.1990 0.030 0.001 0.000
> round(100*固有値/sum(固有値),2)
[1] 86.56 13.07 0.37 0.00

  第2固有値までの累積寄与率は99.63%で、非常に高い。このような場合は第1,2固有値に対応する得点のみを分析すればよい。
  対応分析では、主成分分析や因子分析と同じく、寄与率が高い首位の固有値に対応する行・列の得点の値の大小とそれらの相対関係について分析する。そのため、行の得点と列の得点を同じの画面に配置した散布図 (biplot) がよく用いられている。
  関数 corresp の結果は、関数 plot で第1、2軸の得点の散布図 (biplot) を作成することができる。次にデータセット caith の第1、2軸の得点の散布図を作成するコマンドおよびその結果を図1に示す。

> biplot(caith.ca)

corresp(caith,nf=2)のbiplot

図1 corresp(caith,nf=2) の biplot

  図1から分かるように、目の色が dark の人は髪が black の人が多く、髪色が fair の人は目の色が blue か light の人が多いことが読み取られる。
  関数 biplot に引数 type="row" (あるいは type="col" )を用いることで、行(あるいは列)を基準とした散布図を作成することもできる。
  関数 corresp で計算された行の得点は $rscore、列の得点は $cscore にマトリックス形式で記録されている。

2.2 多重対応分析

(1) 関数 corresp を用いる場合

  説明の便利のため、3つの項目A、B、Cについて10人のアンケート回答の結果を考えよう。項目Aには2つの選択肢(A1、A2)、項目Bには3つの選択肢(B1、B2、B3)、項目Cには3つの選択肢(C1、C2、C3)があるとする。
  例えば、10人の回答結果は表4のようになったとする。表4(a)では、回答の番号を文字列で記入し、表4(b)では回答番号を数字に置き換えて記入している。また、表4のデータは、表5のように項目を選択している場合は「1」、選択していない場合は「0」で記入することも可能である。

表4 10人のアンケート回答の結果 (a)
   項目 
A B C
1 A1 B1 C2
2 A1 B2 C3
3 A1 B3 C1
4 A1 B1 C3
5 A1 B1 C2
6 A2 B2 C1
7 A2 B2 C3
8 A2 B1 C1
9 A2 B3 C1
10 A2 B2 C2

表4 10人のアンケート回答の結果 (b)
  項目
A B C
1 1 1 2
2 1 2 3
3 1 3 1
4 1 1 3
5 1 1 2
6 2 2 1
7 2 2 3
8 2 1 1
9 2 3 1
10 2 2 2

表5 10人のアンケート回答の (0,1) データ
項目A 項目B 項目C
A-1 A-2 B-1 B-2 B-3 C-1 C-2 C-3
1 1 0 1 0 0 0 1 0
2 1 0 0 1 0 0 0 1
3 1 0 0 0 1 1 0 0
4 1 0 1 0 0 0 0 1
5 1 0 1 0 0 0 1 0
6 0 1 0 1 0 1 0 0
7 0 1 0 1 0 0 0 1
8 0 1 1 0 0 1 0 0
9 0 1 0 0 1 1 0 0
10 0 1 0 1 0 0 1 0

  このような表4および表5のデータ形式の対応分析を多重対応分析と呼ぶ。表 4(b) および表5は数量化V類で扱っているデータ形式である。関数 corresp を用いて、多重対応分析を行うためには、データを表5のように「0」「1」形式にすることが必要である。
  関数 corresp を用いた多重対応分析を説明するため、まず表5のデータセットを作成することにする。

> hyou5<-matrix(c(1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,1,1,0,0,1,0,0,0,1,0,0,0,1,1,0,0,1,0,0,1,0,0,0,0,0,1,0,0,0,1,0,0,1,0,1,1,0, 1,0,0,0,1,0,0,0,0,1,0,1,0,1,0,0,1,0,0,0),10,8)
> colnames(hyou5)<-c("項目1.A1","項目1.A2","項目2.B1","項目2.B2","項目2.B3","項目3.C1","項目3.C2","項目3.C3")

  作成したデータ hyou5 の対応分析のコマンドを次に示す。

> (hyou5.ca<-corresp(hyou5,nf=3))

First canonical correlation(s): 0.7740774 0.6915137 0.5773503
Row scores:

[,1] [,2] [,3]
[1,]
-1.3221791 -0.907832620 -0.6800204
[2,]
-0.5536661 1.163884128 1.3600408
[3,]
0.9172655 -1.381667982 1.1900357
[4,]
-1.0956636 -0.008541258 1.3600408
[5,]
-1.3221791 -0.907832620 -0.6800204
[6,]
1.0405163 0.729755881 -0.6800204
[7,]
0.1976570 1.733249149 0.3400102
[8,]
0.4985187 -0.442669504 -0.6800204
[9,]
1.6685886 -0.812302961 0.1700051
[10,]
-0.0288584 0.833957787 -1.7000510

Column scores:

[,1] [,2] [,3]
項目1-A1
-0.8723733 -0.5905856 8.833724e-01
項目1-A2
0.8723733 0.5905856 -8.833724e-01
項目2-B1
-1.0468924 -0.8195340 -2.944575e-01
項目2-B2
0.2117517 1.6127109 -2.944575e-01
項目2-B3
1.6702814 -1.5863538 1.177830e+00
項目3-C1
1.3321954 -0.6893878 8.741354e-16
項目3-C2
-1.1511410 -0.4732166 -1.766745e+00
項目3-C3
-0.6251195 1.3924004 1.766745e+00

  誌面の都合により、biplot(hyou5.ca) 散布図のみを図2 (第1、2得点の散布図) に示す。

corresp(hyou5)のbiplot

図2 corresp(hyou5) の biplot

  表5のようなデータセットを作成するのは若干面倒である。特に質問項目の中に選択肢が多い場合は、列の数が膨大になり、入力する際に間違いやすい。よって、表4のように集計するのが多い。そこで群馬大学社会情報学部青木繁伸氏は、表 4(b) のようなデータを表5のような「0」「1」形式に変換する関数 make.dummy を作成し、次のサイトで公開している。また、そのページではR用の数量化V類のプログラムも公開している。ホームページから数量化V類と関数 corresp の結果が一致することを確認することができる。

  http://aoki2.si.gunma-u.ac.jp/R/qt3.html

(2) 関数 mca を用いる場合

  パッケージ MASS には、多重対応分析のための関数 mca がある。関数 mca の書き式を次に示す。

mca(df, nf ,)

  引数 df はデータフレーム、nf は軸の数で、デフォルトは2になっている。関数 mca は表 4(a) の形式を対象としている。
  R上でデータフレーム形式のデータセットの作成は幾つかの方法がある。そのひとつは、まず空のデータフレームを作成し、データ作成エディタで各セルの値を入力する方法である。
  表 4(a) のデータを hyou4 という名前のデータフレームの作成過程を次に示す。

> hyou4<-data.frame()
> fix(hyou4)

  コマンド fix(hyou4) を実行すると空のデータシートが開かれる。各セルにデータの入力が終ったらシートを閉じることで入力したデータが hyou4 に格納される。
  関数 mca によるデータセット hyou4 の多重対応分析のコマンドおよびその結果を次に示す。

> (hyou4.mca <- mca(hyou4))

Call:
mca(df = hyou4)
Multiple correspondence analysis of 10 cases of 3 factors
Correlations 0.774 0.692 cumulative % explained 38.70 73.28

  コマンド summary(hyou4.mca) で hyou4.mca に格納されている結果のリストが確認できる。正準相関は $d、行の得点は $rs、列の得点は $cs に記録されている。

> hyou4.mca$rs

1 2
1 -0.107883095 -0.0661736894
2 -0.045176338 0.0848377829
3 0.074844205 -0.1007124726
4 -0.089400586 -0.0006225889
5 -0.107883095 -0.0661736894
6 0.084900845 0.0531933287
7 0.016127811 0.1263399093
8 0.040676597 -0.0322670431
9 0.136148354 -0.0592103463
10 -0.002354699 0.0607888088

  関数 corresp の計算結果と関数 mca の結果の各々数値は異なるが、散布図の結果は基本的には同じであることが分かる。関数 mca の得点の biplot を図3に示す。

> biplot(hyou4.mca$rs,hyou4.mca$cs, var.axes = FALSE)

mca(hyou4)のbiplot

図3  mca(hyou4) の biplot

3.その他の対応分析パッケージ

  対応分析には多くのアルゴリズムが提案されている。MASS以外に、パッケージ ade4 や vegan などに対応分析の関数がある。パッケージ ade4 の中に同梱されている対応分析の関数を表6に示す。

表6 パッケージ ade4 における対応分析の関数
関数名 名称
 dudi.coa
 対応分析
 dudi.acm
 多重対応分析
 dudi.fac
 ファジイ対応分析
 dudi.nsc
 非対称(non symmetric)対応分析
 dudi.dec
 偏(decentered)対応分析
 foucart
 K-tableの対応分析

  関数 dudi.coa は関数 corresp、関数 dudi.acm は関数 mca に対応する。計算結果も基本的には同じである。ただし、関数 dudi.coa と関数 dudi.acm では正準相関の代わりに固有値を返す ($eig に格納されている)。得点の散布図は関数 scatter を用いて作成することができる。次にそのコマンド例を示す。関数 dudi.coa では、データフレーム形式しか扱わないことに注意が必要である。

> hyou5<-as.data.frame(hyou5)
> hyou5.coa<-dudi.coa(hyou5,scannf=FALSE,nf=3)
> scatter(hyou5.coa,posieig = "none")

図4 dudi.coa(hyou5)の第1,2軸の散布図

図4 dudi.coa (hyou5) の第1,2軸の散布図

  散布図の軸は次のように指定する。

> scatter(hyou5.coa,xax = 1, yax = 3,posieig = "none")#第1軸と3軸の散布図の作成

  パッケージ vegan には関数 cca、decorana がある。前者 cca は正準対応分析 (Canonical Correspondence Analysis) で、後者 decorana は DCA (Detrended Correspondence Analysis) である。
  対応分析では、図2〜4のように点の散布形態が馬蹄のような形になるケースが多い。この問題を馬蹄問題あるいはアーチ (arch) 問題と呼ぶ。DCAは、アーチ効果を除去する対応分析の方法である。
  また、パッケージ made4 には重み付きの対応分析の関数 dudi.rwcoa (Row weighted Correspondence Analysis) がある。
  また、パッケージ made4 には重み付きの対応分析の関数 dudi.rwcoa (Row weighted Correspondence Analysis) がある。

> library(ade4);library(vegan);
> caith.dca<-decorana(caith)
> s.label(caith.dca$cproj, clab = 1.3)
> s.arrow(caith.dca$rproj, add.pl = TRUE)

図5 decorana(caith)の散布図

図5 decorana (caith) の散布図

参考文献:
[1] 高根芳雄(1995):制約つき主成分分分析法、朝倉書店
[2] 柳井晴夫(1994):多変量データ解析法―理論と応用―、朝倉書店