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

Rと対応分析

1.対応分析とは

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

表1 データ行列Fr×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あるいはPr×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=ZtZ、行の効果はQt=ZZtの固有ベクトル(U,V)をそれぞれ1/√p●j1/√pi●で基準化した値(Dc-1/2UDr-1/2V)を用いる。ここのDcp●jを要素とした対角行列、Drpi●は を要素とした対角行列である。
 Q=ZtZの固有値がλ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

caithの結果

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

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

(caith.ca<-corresp(caith,nf=4))の結果

 ここでは軸数の引数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は数量化Ⅲ類で扱っているデータ形式である。関数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))

 (hyou5.ca<-corresp(hyou5,nf=3))の結果①

 (hyou5.ca<-corresp(hyou5,nf=3))の結果②

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


corresp(hyou5)のbiplot

図2  corresp(hyou5)のbiplot

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

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

> hyou4.mca$rs
 (hyou5.ca<-corresp(hyou5,nf=3))の結果

 関数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):多変量データ解析法―理論と応用―、朝倉書店