統計学といくつかのよしなしごと

つれづれなるままに、日くらし統計学にむかひて、心にうつりゆくよしなし事を、そこはかとなく書きつくれば、あやしうこそものぐるほしけれ。

コレポンは何をやっているのか①数理の概要

告解

 ここ数か月とても忙しく仕事をしていて*1、前の投稿から半年以上経ってしまった。とはいえ何よりも自分の知識や経験の定着を一番の目的にブログを書いているので、これだけ空いたのは怠慢でしかないんだよなぁ。。もう少し頑張ろうと思います。

マーケティングリサーチでよくやるいわゆるコレポン

 気を取り直して。マーケティングリサーチではコレスポンデンス分析*2、いわゆるコレポンがよく用いられる。マーケティングはもちろん、戦略や戦術を考えるといった際に僕たちは二軸に切ったマップを頻繁に利用する。会社であれ製品・サービスであれ顧客ニーズであれ、切れ味の良さそうな軸を定義してその上にそれらを置いてやれば、何となく洞察が得られた気になる。実際関心のある状況を簡略化、可視化することは思考の大きな助けになるし、その点マッピングはとても有用である。しかし経験的直感だけでマッピングするのは心許ないから、データに基づいて確からしマッピングしたい。これに応える手法の一つがコレスポンデンス分析である。
 例として、あるモデル事務所が目新しいモデルを打ち出したいと思っており、目新しさとは何か検討するために従来のモデルの外見的特徴をマッピングする、という架空のストーリーを考える。外見の9割は目と髪の色で決まると考え、各モデルにこれらの値を付与したデータセットを作成した*3。これに対しコレスポンデンス分析は以下のようなマップを出力する。図の黒文字が目の色、赤文字が髪の色である。

> 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
> # 表側:目の色
> # 表頭:髪の色
> # 各セルの値は該当する人数

f:id:yanbow221:20180716230923j:plain

 これを見ると目の色と髪の色の組み合わせとして[dark , black][blue or light , fair]などのモデルが多いことが分かる。プロットされた点から何となくモデルの外見的特徴は「暗め」「明るめ」「標準(medium)」の3タイプに分かれそうな気がするし、少なくとも横軸については「明るい-暗い」といった意味付けをしたくなる。このようにコレスポンデンス分析は関心のあるデータ項目を数理的にプロットしてくれるため、ビジネスを検討する上での重要かつ使い勝手の良いツールとなるのである。
 しかし、そもそもこの横軸縦軸は何を表しているのか、なぜその2軸が選ばれたのか、なぜ異なる意味合いを持つ項目を同じ図に布置できるのか、軸に意味付けをすることは正当か、正当ならば意味付けに際して注意するべきことはあるか、などなどの疑問が湧く。分かりやすい可視化で容易に利用できる手法であるからこそ、分析者は手法に対する理解を深め、その後のビジネス上の検討を誤った方向にリードしてしまわないように注意する必要があるよなぁ。と思って以下の本で一通り勉強しなおしたので備忘としてブログに残しておくことにする。

対応分析入門 原理から応用まで 解説◆Rで検算しながら理解する

対応分析入門 原理から応用まで 解説◆Rで検算しながら理解する

 Rでの検算の章もあって理解の助けになるが、文章に癖があるためかややとっつきにくいので以下リンクの「対応分析法・数量化Ⅲ類の考え方」等で補完するのが良いかも。
http://www.wordminer.org/wp-content/uploads/2013/04/63_18.pdf

コレスポンデンス分析の数理の概要

 コレスポンデンス分析は対応分析とも言われるように、先の例で見たように表側の項目と表頭の項目を対応させ一つの図に同時に布置する分析手法である。イメージはこんな感じ。
f:id:yanbow221:20180717102910j:plain

  1. 行、列それぞれの方向に対してデータを変換する
  2. 変換したデータについて分析してそれぞれプロットする
  3. 重ね合わせる

 3ステップで直感的にはとても分かりやすい。では具体的にどのようにこれを行っているのか実際に計算しながら見ていくことにする。

データセット

 上でリンクしたPDF「対応分析法・数量化Ⅲ類の考え方」に記載されている以下の架空データを使う。10サンプルについてある商品の「好きな銘柄」を選んでもらったという設定。なので「好き=1, 好きでない=0」の二値データになっている。

## データセット
F <- matrix(c(0, 0, 0, 0, 0, 1, 1, 0, 1, 1,
              1, 0, 0, 1, 1, 1, 1, 0, 1, 0,
              0, 0, 1, 1, 1, 1, 0, 1, 0, 0,
              0, 0, 0, 0, 0, 0, 1, 0, 0, 1,
              1, 0, 0, 1, 0, 1, 1, 0, 1, 1,
              1, 1, 1, 1, 1, 0, 0, 1, 0, 0),
            nrow = 10, ncol = 6, byrow = F)
dimnames(F) <- list(サンプル = c("サンプル1", "サンプル2", "サンプル3", "サンプル4", "サンプル5",
                             "サンプル6", "サンプル7", "サンプル8", "サンプル9", "サンプル10"),
                        銘柄 = c("銘柄A", "銘柄B", "銘柄C", "銘柄D", "銘柄E", "銘柄F"))

> F
            銘柄
サンプル     銘柄A 銘柄B 銘柄C 銘柄D 銘柄E 銘柄F
  サンプル1      0     1     0     0     1     1
  サンプル2      0     0     0     0     0     1
  サンプル3      0     0     1     0     0     1
  サンプル4      0     1     1     0     1     1
  サンプル5      0     1     1     0     0     1
  サンプル6      1     1     1     0     1     0
  サンプル7      1     1     0     1     1     0
  サンプル8      0     0     1     0     0     1
  サンプル9      1     1     0     0     1     0
  サンプル10     1     0     0     1     1     0
ステップ1:行、列それぞれの方向に対してデータを変換する

 まずはマッピングするためのデータを行、列それぞれに作っていくステップ。ここでもいくつか踏む手順があるので、順番に見ていく。
1.元のデータセットをの確率行列を作成する
 第1段階として同時確率分布{\displaystyle \bf P_{ij}}、行方向の周辺確率分布を対角要素にした対角行列{\displaystyle \bf P_I}、列方向の周辺確率分布を対角要素にした対角行列{\displaystyle \bf P_J}を用意する。

{\displaystyle 
\begin{eqnarray}
{\bf P_{ij}} &=& \left(p_{ij}\right) = \left(\frac{f_{ij}}{N}\right),\ \ \ N = \sum^m_{i=1}\sum^n_{j=1}f_{ij}
\\
{\bf P_I} &=& diag \left( p_{i+} \right) = \left(\frac{f_{i+}}{N}\right) = \left(\frac{\sum^n_{j=1} f_{ij}}{N}\right)
\\
{\bf P_J} &=& diag \left( p_{+j} \right) = \left(\frac{f_{+j}}{N}\right) = \left(\frac{\sum^m_{i=1} f_{ij}}{N}\right)
\end{eqnarray}
}

2.確率行列を使って行、列それぞれのプロファイルを作成する
 ここで行、列それぞれに対応する行列が作られる。要は行方向に見たときの条件付き確率分布、列方向に見たときの条件付き確率分布のことだが、対応分析においてはそれぞれ行プロファイル列プロファイルと呼ばれる。定式化するとこんな感じ。

{\displaystyle 
\begin{eqnarray}
{\bf N_I} &=& \left( q_{ij} \right) &=& \left( \frac{p_{ij}}{p_{i+}} \right) \ \ \ (行プロファイル)
\\
{\bf N_J} &=& \left( q_{ij}^* \right) &=& \left( \frac{p_{ij}}{p_{+j}} \right) \ \ \ (列プロファイル)
\end{eqnarray}
}

3.ストレッチプロファイルを作成し平均値で中心化する
 上記で作成した行列に対して距離を求めた場合はユークリッド距離となるが、独立性の検定に用いられるカイ二乗統計量に結び付けるためカイ二乗距離になるようデータを変換したものをストレッチプロファイルと呼ぶ(右辺第1項)。こうして作った行列から平均値(右辺第2項)を引いて中心化した式が以下で、これが対応分析に利用する行列になる。

{\displaystyle 
\begin{eqnarray}
{\bf X} &=& \left( x_{ij} \right) &=& \left( \frac{q_{ij}}{\sqrt{p_{+j}}} - \sqrt{p_{+j}} \right)
\\
{\bf X^*} &=& \left( x_{ij}^* \right) &=& \left( \frac{q_{ij}^*}{\sqrt{p_{i+}}} - \sqrt{p_{i+}} \right)
\end{eqnarray}
}

 なお、この式に関する詳細な説明は以下の「対応分析法の基本的な考え方(数理の要点)」が詳しい。
http://wordminer.org/wp-content/uploads/2013/04/2d3ed444f1290ae0e0992a05b64a5817.pdf

 ここまでの作業を書いたコードは以下。

# 確率行列
P <- F / sum(F)

# 行の確率要素(P_iの対角要素)と対角行列化
pi <- matrix(apply(P, 1, sum))
Pi <- diag(pi[, 1])

# 列の確率要素(P_jの対角要素)と対角行列化
pj <- matrix(apply(P, 2, sum))
Pj <- diag(pj[, 1])

# 行プロファイル
Ni <- diag(1 / pi[, 1]) %*% P

# カイ二乗距離を求めるためのストレッチプロファイル
Ni2 <- Ni %*% diag(1 / sqrt(pj[, 1]))
dist(Ni2) # ストレッチしているためdist関数で求める距離はカイ二乗距離になる

# 平均値(列質量の平方根)による中心化
Xi <- sweep(Ni2, 2, sqrt(pj), FUN = "-") # 対応分析における基本のデータ行列①(行)


# 列プロファイル
Nj <- diag(1 / pj[, 1]) %*% t(P)

# カイ二乗距離を求めるためのストレッチプロファイル
Nj2 <- Nj %*% diag(1 / sqrt(pi[, 1]))
dist(Nj2) # ストレッチしているためdist関数で求める距離はカイ二乗距離になる

# 平均値(行質量の平方根)による中心化
Xj <- sweep(Nj2, 2, sqrt(pi), FUN = "-") # 対応分析における基本のデータ行列②(列)
ステップ2:変換したデータについて分析してそれぞれプロットする

 こうして作られた行列{\displaystyle \bf X}{\displaystyle \bf X^*}に対して以下の操作を加えてやるとそれぞれの分散共分散行列になる。

{\displaystyle 
\begin{eqnarray}
{\bf V} &=& {}^t\! {\bf X} {\bf P_I} {\bf X}
\\
{\bf V^*} &=& {}^t\! {\bf X^*} {\bf P_J} {\bf X^*}
\end{eqnarray}
}

 これらの行列についてそれぞれ固有値問題を解けば固有値固有ベクトルが求まるが、これはまさに主成分分析である。つまり対応分析は行方向と列方向のそれぞれに対して主成分分析を行っていることになる。
 ちなみに主成分分析が固有値問題に帰着する理由については毎度のことながら永田先生の本がとても分かりやすい。

多変量解析法入門 (ライブラリ新数学大系)

多変量解析法入門 (ライブラリ新数学大系)

 ブログだと以下が丁寧に書かれてある。
blog.aidemy.net

 こうして求めた固有ベクトルで元のデータ(ストレッチプロファイル)を線形変換してやることで目当ての座標系が得られる。コードはこんな感じ。

### 固有値問題に持ち込む
## 行方向
# 分散共分散行列
Vxi <- t(Xi) %*% Pi %*% Xi

# 固有値問題を解く
result_Xi <- eigen(Vxi)

# 結果の確認
result_Xi$values # 固有値
result_Xi$vectors # 固有ベクトル

# 固有ベクトルによる座標変換
rscore <- Ni2 %*% result_Xi$vectors

## 列方向
# 分散共分散行列
Vxj <- t(Xj) %*% Pj %*% Xj

# 固有値問題を解く
result_Xj <- eigen(Vxj)

# 結果の確認
result_Xj$values # 固有値
result_Xj$vectors # 固有ベクトル

# 固有ベクトルによる座標変換
cscore <- Nj2 %*% result_Xj$vectors

 なお固有値の個数は{\displaystyle min\{m, n\} - 1}で、今回の場合は5個となる。1を引いているのは元々プロファイルが確率行列で和が行和ないし列和が1となる制約条件による。解いた固有値固有ベクトルを確認してみると確かに5成分で表現できていそう。
 ところでこれまでは行と列個別に見てきたが、これら2つには双対性という性質があり、列成分スコアは行成分スコアを用いて、行成分スコアは列成分スコアを用いて求めることができる。算出式を遷移方程式と言って、行成分スコア(rscore)を{\displaystyle \bf Z_{ik}}、列成分スコア(cscore)を{\displaystyle \bf Z_{jk}^*}とすれば(kは固有値/成分に該当)、

{\displaystyle 
\begin{eqnarray}
{\bf Z_{ik}} &=& \frac{1}{\sqrt{λ_k}} \sum^n_{j=1} \left( \frac{p_{ij}}{p_{i+}} \right) z_{jk}^*
\\
{\bf Z_{ik}^*} &=& \frac{1}{\sqrt{λ_k}} \sum^m_{i=1} \left( \frac{p_{ij}}{p_{+j}} \right) z_{ik}^*
\end{eqnarray}
}

で定式化される。コードは以下。とりあえずマッピングに際しては列成分スコアを遷移方程式で求めなおして利用することにする。

# rscoreによる遷移方程式を用いてcscoreを求めなおす
rscore2 <- rscore[, 1:5] # 対象の5成分を抽出
C <- diag(1 / pj[, 1]) %*% t(P) # 列プロファイルの再算出
Dm1 <- diag(1 / sqrt(result_Xi$values[1:5]))
cscore2 <- C %*% rscore2 %*% Dm1 # 遷移方程式

# 成分スコアの確認
dimnames(rscore2) <- list(サンプル = c("サンプル1", "サンプル2", "サンプル3", "サンプル4", "サンプル5",
                                   "サンプル6", "サンプル7", "サンプル8", "サンプル9", "サンプル10"),
                              成分スコア = c("第1成分", "第2成分", "第3成分", "第4成分", "第5成分"))

dimnames(cscore2) <- list(銘柄 = c("銘柄A", "銘柄B", "銘柄C", "銘柄D", "銘柄E", "銘柄F"),
                            成分スコア = c("第1成分", "第2成分", "第3成分", "第4成分", "第5成分"))
rscore2 # 行成分スコア
cscore2 # 列成分スコア

 

ステップ3:重ね合わせる

 ひとまず横軸に第1成分、縦軸に第2成分をとって可視化してみると、こんな感じ。

# 第1成分、第2成分による同時布置図
library(ggplot2)

dat_map <- rbind(rscore2, cscore2)
dat_map <- transform(dat_map, label = c(rep("サンプル", 10), rep("銘柄", 6)))

map1 <- ggplot(dat_map, aes(x =1成分, y =2成分, color = label)) + 
  xlim(-1.5, 1.5) + ylim(-1.5, 1.5) +
  geom_hline(yintercept = 0, color = "black") +
  geom_vline(xintercept = 0, color = "black") +
  geom_point(size = 5) +
  scale_color_manual(
    values = c("サンプル" = "tomato", "銘柄" = "skyblue")
  ) +
  geom_text(aes(label = rownames(dat_map)), vjust = -1, size = 5) +
  theme_bw()

map1

f:id:yanbow221:20180717203905j:plain

これで対応分析の主要アウトプットである同時布置図を描くことができた。

釈然としない

 対応分析の流れを追ってみたが、行成分スコアと列成分スコアが同じ次元にマッピングされる正当性について明瞭に示されている感じはしない。実は特殊なケースを除いてこれら2つのスコアが同次元にはならず、かい離があることが指摘されていて、同じ空間にマッピングされるのは解釈に便利だから、ということらしい。さらに、実際計算してみると分かるが、上で紹介した「対応分析法・数量化Ⅲ類の考え方」に記載されている行成分スコア・列成分スコアとRで計算したスコアは一部符号が反転している成分がある。これは固有ベクトルの持つ性質*4とそれに対する利用ソフトウェアのアルゴリズムが異なるためで、つまり同じデータセットに対して対応分析を行ってもソフトウェアによってマップが視覚的に異なる可能性があることを示している。次元が異なるならうかつに軸に意味付けすることができないし、やはり対応分析の結果の解釈については注意が必要であることを裏付けている。

 長くなってしまったで対応分析の結果をどのように評価していけば良いかについては次のブログでまとめることにします。

*1:お客さんとの雑談の中でPJメンバーの一人がリアルに「欲しがりません勝つまでは」と言うくらい大変だった

*2:対応分析とも言われる

*3:Rのパッケージ「MASS」で提供されているデータセット「caith」を使うためこのようなあるかどうかも分からない設定にした。実際にはcaithはイギリス在住者5,387名に対する調査結果らしい

*4:固有ベクトルvに対してスカラーkを掛けたkvもまた固有ベクトルであるという性質。このkをどう定めるかがソフトウェアによって異なるだろうため、符号の反転が起こる