爆速習・コレスポンデンス分析（対応分析）

<details><summary>Disclaimer</summary>
・この記事は筆者の所属する組織・団体とは一切関係がありません。<br>
・この記事の内容を利用することによって生じた損失等について、筆者は一切責任を負いません。
</details>

・この記事ではコレスポンデンス分析（対応分析）について学べます。
・線形代数の基礎的な知識がある人向けの説明が含まれます。

## コレスポンデンス分析を爆速で学ぶ
コレスポンデンス分析はマトリクス型のカウントデータに対して適応可能な手法です。
順番に説明していきましょう。

#### 架空の例
70人の男性にクワガタ虫に対する意識調査を行いました。
具体的には、五種類のクワガタ「ノコギリクワガタ・ミヤマクワガタ・コクワガタ・オオクワガタ・ヒラタクワガタ」それぞれに対して、「かっこいい・かわいい・強そう・弱そう」それぞれに当てはまるかを聞きました。

例えば、ある男性は
ノコギリ：かっこいい
ミヤマ　：かっこいい・強そう
コクワ　：かわいい・弱そう
オオクワ：かっこいい・かわいい・強そう
ヒラタ　：かっこいい・強そう
のように回答しています。こういった回答が70人分あります。

これらを集計した結果が以下の表になります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3784717/babef006-6c62-4174-1c71-81d1efd3fce1.png)

やはり、オオクワ・ヒラタは「かっこいい」イメージですね。コクワは「かわいい」イメージが強いみたいです。
「かっこいい・強そう」はプラスに相関していそうですね。「かっこいい」クワガタは「強そう」ということでしょう。
「強そう・弱そう」はマイナスに相関していそうですね。「強そう」なクワガタは「弱そう」ではない、という当たり前の結果です。

さて、**それぞれのクワガタはどのような特徴を持っているでしょうか？　各クワガタのポジションを調べたい！**
……そのための分析手法がコレスポンデンス分析となります。

今回の例はクワガタのイメージ調査ですが、例えば「企業A,B,C……のブランドイメージ」や、「製品A,B,C……のユーザー評価」でも利用できる手法となります。

#### コレスポンデンス分析
各クワガタのポジションを調べたいわけですが、**行項目（クワガタ）に数値を割り振る**ことを考えてみましょう。
例えば「ノコギリ・ミヤマ・コクワ・オオクワ・ヒラタ」に対して、それぞれ「-5,-4,0,5,3」のように数値を割り振れば、「ノコギリとミヤマは似ている」「オオクワとヒラタは似ている」「ノコギリ→ミヤマ→コクワ→ヒラタ→オオクワ、という順序がある」ということになりますよね。
当然ながら数値は適当に割り振っても意味がないので、データに基づいて行いたいです。
**数値の割り振りを、列項目（各クワガタのもつイメージ）を基準にして行いたい**わけです！

例えば列項目の一つ「かっこよさ」をそのまま引っ張ってきて割り振るなら、
「ノコギリ・ミヤマ・コクワ・オオクワ・ヒラタ」はそれぞれ「35・32・23・60・55」となりますね。「コクワ→ミヤマ→ノコギリ→ヒラタ→オオクワ」という順序が各クワガタに存在していることがわかります。ここでの順序はもちろん、「かっこよさの順序」ですね。

このままでは列項目の一つ「かっこよさ」しか使っていません。列項目は全部で四つあるのですから、列項目を全て考慮した数値を行項目に割り振りたいですね。

さて、**列項目にも数値を割り振る**ことを考えてみましょう。
行項目（クワガタ）に割り振る数値を$x$とします。列項目（イメージ）に割り振る数値を$y$とします。
具体的には、「ノコギリ・ミヤマ・コクワ・オオクワ・ヒラタ」それぞれに対して割り振る数値をそれぞれ $x_1,x_2,x_3,x_4,x_5$ として、「かっこいい・かわいい・強そう・弱そう」それぞれに対して割り振る数値を $y_1,y_2,y_3,y_4$ とします。
そうすると、**このマトリクスデータは$(x_i,y_j)$というデータがたくさん観測されたものを集計した結果である、と見ることができます**。もう少し詳しく説明しましょう。

以下がマトリクスの再掲です。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3784717/babef006-6c62-4174-1c71-81d1efd3fce1.png)
「ノコギリ」「かっこいい」のセルは35ですから、35人が「ノコギリ」は「かっこいい」と評価したことになります。**これを$(x_1,y_1)$が35回観測された、と考えてみます**。
このマトリクスの全てのセルを足し上げると640になるので、**$(x,y)$というデータが全部で640個観測されていて、そのうち35個は$(x,y)=(x_1,y_1)$だった、と考える**わけですね。

さて、具体的にどうやって$(x_i,y_j)$の値を決めていけば良いでしょうか。
いま、$(x,y)$という二変量のデータが640件観測されています。
ということは、**サンプルサイズ640で相関係数を計算することができます**よね！
(※もちろん、まだ具体的な数値は出せません。具体的な$x,y$の数値を割り振ってないので)

**★この相関係数を最大化するように$x,y$の値を決めてみます★**
これがコレスポンデンス分析の核となる考え方です。
相関係数が大きくなるように$x,y$の値を割り振れば、「$x$が大きくなるほど$y$が大きくなる」ような関係性をつくることができます。$x$は行項目、$y$は列項目に割り振る値でしたから、これはすなわち行項目と列項目を密接に結びつけるような数値を、行・列それぞれに割り振ることを意味しています。
このようにして決めた$x$を使って各クワガタをグルーピング・順序化してやれば、($x$と$y$は強く相関するのだから)それは自動的に列項目に基づいたグルーピング・順序化になっているわけです。つまり、各クワガタの持つイメージ（列項目）に基づいて各クワガタに数値を割り振ることに成功しているわけです！

#### 数理的な説明
とりあえず$x_{vec} = (x_1,\cdots,x_5)^\prime , y_{vec}=(y_1,\cdots,y_4)^\prime$として縦ベクトルを定義しておきましょう。${}^\prime$は転置記号です。

マトリクスを$N_{all}$で表しましょう。今回の例であれば
$$
N_{all} = \left( \begin{array}{cccc} 
35	&15	&31	&27 \\\ 
32	&13&	38	&22\\\ 
23	&51&	9	&49\\\ 
60&	22&	65&	12\\\ 
55&	30&	46&	5
\end{array} \right)
$$
となります。
これの行和を対角に並べた対角行列を$N_{row}$とし、列和を対角に並べた対角行列を$N_{col}$とします。また、行和を並べた縦ベクトルを$n_{row}$とし、列和を並べた縦ベクトルを$n_{col}$とします。すると、$x$の総和は$n_{row}^\prime x_{vec}$で表せますね。同様に$y$の総和は$n_{col}^\prime y_{vec}$で表せます。

さて、このまま相関係数最適化をしようとしても最適な組み合わせが無限個だけ出てきてしまうので、$x,y$がとれる値の範囲に制約を付けます。
具体的には「総和がゼロ（したがって平均もゼロ）」という制約を付けます。つまり、
$$
n_{row}^\prime x_{vec}=0,  \ \ \ n_{col}^\prime y_{vec} = 0
$$
とします。

さて、$x,y$とも平均がゼロなわけですから、相関係数を$r$とすると、
$$
r = \frac{y_{vec}^\prime N_{all} x_{vec} }{\sqrt{x_{vec}^\prime N_{row} x_{vec} \ \ y_{vec}^\prime N_{col} y_{vec} } }
$$
となりますね。これを制約付き最適化するわけですから、ラグランジュの未定乗数法より最終的に(途中式は省略)、
$$
N_{all} N_{col}^{-1}N_{all}^\prime x_{vec} = r^2 N_{row} x_{vec}
$$
の解を求めればよいことがわかり、ここで$w_{vec} = N_{row}^{1/2} x_{vec}$とすれば
$$
N_{row}^{-1/2} N_{all} N_{col}^{-1}N_{all}^\prime N_{row}^{-1/2} w_{vec} = r^2 w_{vec}
$$
の解を求めればよいとわかります。左辺の行列がゴチャゴチャしてて分かりにくいので単に$\Lambda = N_{row}^{-1/2} N_{all} N_{col}^{-1}N_{all}^\prime N_{row}^{-1/2}$とまとめてしまえば、
$$
\Lambda w_{vec} = r^2 w_{vec}
$$
の解を求めればよく、これは単純に$\Lambda$の固有値問題です。

ここで、第一固有値(一番大きい固有値)$r^2 = 1$ および対応する固有ベクトルは自明な解なので考えないことにして、第二固有値に対応する固有ベクトル$w^{(2)}_{vec}$を求めてやり、それを更に$x^{(2)} _{vec} = N^{-1/2} _{row} w^{(2)} _{vec}$で$x$のスケールに戻すことにします。

**$x^{(2)} _{vec}$が、求めたかった「行項目に割り振る数値」となります**。
また、第三固有ベクトル$w^{(3)}_{vec}$に対応する$x^{(3)} _{vec}$も考えることができますね。（これは相関係数最適化問題の局所解です）

### 実装結果
上述のマトリクスを自分で作って"ダミーデータ.csv"ファイルに保存してください。
そうすると、以下のコードでコレスポンデンス分析が実行できます。

```R:
Data <- read.csv(file = "ダミーデータ.csv", row.names = 1)
N_all <- as.matrix(Data)
N_row <- diag(rowSums(N_all))
N_col <- diag(colSums(N_all))
Lambda.row <- 
  solve(N_row)^(1/2) %*% N_all %*% solve(N_col) %*% t(N_all) %*% solve(N_row)^(1/2)
res.eigen <- eigen(Lambda.row)
eigen.values.row <- res.eigen$values
eigen.vecs.row <- res.eigen$vectors
plot(eigen.vecs.row[,2], rep(0,5) , pch = 16 , cex = 1.5)
abline(h = 0 , lty = 2)
text(x = eigen.vecs.row[,2] , y = 0.2 , rownames(Data))
```
$x^{(2)} _{vec}$をプロットした結果が以下の図になります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3784717/fbb47d0e-ef0d-b40c-e059-3b4f14a131e3.png)
これを見ると、コクワとオオクワは対極の位置にあり、「コクワ→ノコギリ→ミヤマ→ヒラタ→オオクワ」という「順序」を見出すことができますね。この数値は「弱そう」から「強そう」への変化、「かわいい」から「かっこいい」への変化を表していると考えられます。

さて、第三固有ベクトルも考えてみましょう。
```R:
plot(eigen.vecs.row[,2], eigen.vecs.row[,3] , pch = 16 , cex = 1.5 , 
     ylim = c(min(eigen.vecs.row[,3])-0.2,max(eigen.vecs.row[,3])+0.2))
abline(v = 0, h = 0 , lty = 2)
text(x = eigen.vecs.row[,2] , y = eigen.vecs.row[,3]+0.1 , rownames(Data))
```

$x^{(2)} _{vec}$を横軸、$x^{(3)} _{vec}$を縦軸にプロットした結果が以下の図になります。
**※これがいわゆる「ポジショニングマップ」になります**。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3784717/2bd0d30d-9366-0dbd-3dca-40c6bbfeb0a5.png)
少し面白い結果となりました。
第一指標(横軸)では正反対だったコクワとオオクワが、第二指標(縦軸)では「似ている」と判断されています。実際、コクワガタとオオクワガタは同じドルクス族に属するクワガタであり、姿かたちはとてもよく似ています(大きさ以外は！）。そういった点で、オオクワガタとコクワガタが持つイメージは似ている側面もあるということでしょう。
ノコギリとミヤマは横軸でも縦軸でも非常に近い位置にあります。この二種はかなり似たポジションにあるとわかります。

## 参考文献
『分割表の統計解析 二元表から多元表まで』 宮川雅巳・青木敏, 朝倉書店, 2018