対応分析の結果の解釈について

はじめまして、データ分析部の青木です。

今日は多変量解析法の一つである対応分析(コレスポンデンス分析)を題材として、分析結果の解釈の問題についてお話ししたいと思います。

対応分析はクロス集計表の行と列のカテゴリーの関係性を視覚的に把握するための分析手法です。マーケティング・リサーチの分野でよく用いられるため、皆さんの中にも馴染みのある方が多いかと思います。

まずはRでの対応分析の実行方法を簡単に紹介したいと思います。Rで対応分析といえばMASSパッケージのcorresp関数がポピュラーかと思いますが、今回はcaパッケージを利用したいと思います。

install.packages("ca")
library("ca")

サンプルデータとしてcaパッケージに含まれているauthorデータを使うことにします。authorデータは、12種類の小説に対する本文の一部(固有名詞を除く)からアルファベットの出現頻度をカウントしmatrix形式でまとめたものです。

data(author)
author

以下のデータが表示されると思います。(合計の列は私が付けたものです。)
author1
ここでは小説ごとにアルファベットの出現頻度(相対度数)に差があるかどうかを調べたいとします。カイ二乗適合度検定を行なうのも一つの手ですが、これだけサンプルサイズがあるとたいていの場合は棄却されてしまうのであまり意味がありませんし、そもそも「どのような差があるのか?」を把握することができません。以下の相対度数の表をじっと眺めてどのような差があるかを調べるのもたいへんです。

author / rowSums(author)

author_row_profile

そこで対応分析の出番となります。使い方は簡単です。(下のコードの1行目と2行目で後々の見易さのために行カテゴリーの名前から題名を取り除いていますが、もちろんこの処理は必須ではありません。)

author.names <- rownames(author)
rownames(author) <- sub(".+\\((\\w+)\\)", "\\1", author.names)
author.ca <- ca(author)
plot(author.ca, map="rowprincipal")

今回は行カテゴリーの相対度数に興味があるのでplot関数の引数でmap=”rowprincipal”と指定していますが、実際はデフォルト値のmap=”symmetric”を設定することが多いかと思います。以下のような図がプロットされると思います。
row_principal_2
列カテゴリーに対して行カテゴリーがほぼ原点に集中して潰れてしまっています。これは小説間でのアルファベットの出現頻度の偏りが小さいことを示しています。が、このままだと相対的な違いが見づらいためmap=”symmetric”(デフォルト値なので指定しなくてもよい)として再プロットします。

plot(author.ca, map="symmetric")

symmetric_2
上図は平面上で近くにプロットされたカテゴリー同士の関係性が強いことを意味しています。どの小説とどのアルファベットの関係性が強いのかも気になりますが、それよりもなんとなく同じ著者の小説同士が近くにプロットされているような気がします。さらに詳細に見るために行カテゴリーのみをプロットすることにします。

plot(author.ca, what=c("all", "none"))

row_only_2
ggplot2パッケージを使うとより分かり易くプロットできます。

install.packages("ggplot2")
library("ggplot2")
rowcoord <- author.ca$rowcoord %*% diag(author.ca$sv)
colnames(rowcoord) <- colnames(author.ca$rowcoord)
rowcoord.df <- data.frame(rowcoord,
                          Author=rownames(author),
                          row.names=1:12)
ggplot(data=rowcoord.df, aes(x=Dim1, y=Dim2, label=Author)) +
  coord_fixed() +
  geom_hline(yintercept = 0, colour = "gray70") +
  geom_vline(xintercept = 0, colour = "gray70") +
  geom_text(aes(colour=Author))

row_only_5
やはり同じ著者の小説同士は近くにプロットされているようです。つまり、著者ごとにアルファベットの出現頻度に偏りがありそうです。この結果は偶然なのでしょうか?例えば以下の文献ではpermutation test(並べ替え検定)を行なうことで、この結果の有意性を示しています。
Correspondence Analysis in Practice, Second Editionの10章と25章

 

では次に、同じcaパッケージの中にあるwg93データに対して対応分析を適用してみたいと思います。

data(wg93)
head(wg93)
  A B C D sex age edu
1 2 3 4 3   2   2   3
2 3 4 2 3   1   3   4
3 2 3 2 4   2   3   2
4 2 2 2 2   1   2   3
5 3 3 3 3   1   5   2
6 3 4 4 5   1   3   2

このデータはInternational Social Survey Program on Environment 1993という調査データの一部で、科学に対する意識調査のための質問4つ(5段階評価で回答)と性別・年代・教育といったデモグラ属性からなるデータです。ここでは対応分析を適用するために上記データの一部(質問AとB)だけを使うことにします。分割表の形にまとめます。

wg93.table <- with(wg93, table(A, B))
wg93.table
   B
A    1  2  3  4  5
  1 27 28 30 22 12
  2 38 74 84 96 30
  3  3 48 63 73 17
  4  3 21 23 79 52
  5  0  3  5 11 29

前の例と同様にca関数で対応分析を行ないます。

wg93.table.ca <- ca(wg93.table)
plot(wg93.table.ca)

wg_2
上図のような結果が得られました。行と列それぞれ同じ番号のカテゴリーが近くにプロットされていることから、この二つの質問には正の相関があることが分かります。ところで、この図の横軸(第1次元)と縦軸(第2次元)はどのように解釈したらよいでしょうか?横軸は1(強くそう思う)~5(全くそう思わない)をそのままストレートに表わしているということでよいかと思いますが、縦軸はどうでしょうか?

実は今回のような順序尺度のデータに対して対応分析を適用すると上図のような放物線状のプロットが得られることがよくあります。この現象はarch effectやhorseshoe effectと呼ばれており、その発生のメカニズムが色々と研究されています。ここでは以下のOtsu (2003)、大津(2007)による考察を辿ってみたいと思います。

まず、下記のように2次元正規分布をもとにした2次元分割表を作成します。

# rmvnorm関数を使うためにmvtnormパッケージを読み込む
install.packages("mvtnorm")
library("mvtnorm")
n <- 10000
mean <- c(0, 0)
sigma <- matrix(c(1.0, 0.7, 0.7, 1.0), nrow=2)
# 2次元正規分布に従う乱数を生成
x <- rmvnorm(n, mean=mean, sigma=sigma)
t <- seq(1, 7, 2) / 3
thresholds <- c(-Inf, -rev(t), t, Inf)
# 連続変数を離散化して分割表にするための閾値
thresholds
 [1]       -Inf -2.3333333 -1.6666667 -1.0000000 -0.3333333
 [6]  0.3333333  1.0000000  1.6666667  2.3333333        Inf
norm.matrix <- matrix(0, 9, 9, dimnames=list(1:9, 1:9))
# 各ビンのサンプルサイズを集計
for (i in 1:9) {
  for (j in 1:9) {
    norm.matrix[i, j] <- sum(thresholds[i] < x[, 1] &
                             x[, 1] <= thresholds[i + 1] &
                             thresholds[j] < x[, 2] &
                             x[, 2] <= thresholds[j + 1])
  }
}
norm.matrix
   1   2   3   4   5   6   7   8  9
1 24  31  29   8   3   1   0   0  0
2 33  93 146  96  25   1   0   0  0
3 22 146 345 370 207  49   3   1  0
4  8  92 386 705 637 227  44   5  0
5  3  17 193 594 979 628 185  20  1
6  0   6  40 229 590 743 335  87  8
7  0   0   2  41 220 388 328 124 23
8  0   0   0   6  35 100 121  94 32
9  0   0   0   1   1   5  27  29 28

そしてこの分割表に対応分析を適用してみましょう。きれいな曲線が得られます。1次元と2次元のプロットでは2次曲線が、1次元と3次元のプロットでは3次曲線が得られています。

norm.matrix.ca <- ca(norm.matrix)
plot(norm.matrix.ca)
plot(norm.matrix.ca, dim=c(1, 3))

horseshoe_12_2
horseshoe_13_2
元のデータ構造が1次元であるにも関わらずこのような多次元の解が得られる理由として、上記文献では正規分布の密度関数とエルミート多項式の関係性に触れています。

H_k(x)k次のエルミート多項式に適当な正規化項をかけたもの、f(x)を標準正規分布の密度関数とすると、{H_k(x)}はf(x)を重みとした場合の正規直交関数系となります。またf(x,y)を相関係数が\rhoである2次元の標準正規分布の確率密度関数とすると、以下の式が成り立ちます。
\int H_i(x)f(x,y)H_j(y)dxdy=\delta_{ij}\rho^i

このことから2次元分割表の背後に2次元正規分布を仮定した場合、高次の多項式が(固有値\rho^{2i}に対応する)解として現れることになります。(厳密には離散化処理があるため近似的にこのような関係性が成立するということになります。)

では逆に、上記のような放物線が得られた場合、もとのデータの構造は1次元であると言えるのでしょうか?答えは残念ながらNoで、そうとは限らないことが興味深い数値例と共に上記文献で紹介されています。

今回は対応分析を例に結果の解釈の問題について見てみました。対応分析に限らずどの手法でも結果の解釈にはとても難しい部分があり、よくよく注意して結果を解釈しないと思わぬ落とし穴にはまってしまうことがあるかと思います。

最後になりましたが、前回の更新でもお知らせしました通り弊社は2015年2月19日に東京証券取引所マザーズ市場へ新規上場をしまして、更なる事業拡大に向けてデータ分析者を募集しております。おかげさまで多数ご応募を頂いておりますがまだまだ積極採用中ですので、ご興味ある方は求人ページをご覧頂ければと思います。

また、私と今井が所属しているデータ分析部R&Dセクションでは、機械学習系の研究室に所属している学生のアルバイトも積極的に募集しております。機械学習を用いた実データの分析に興味がある方は、こちらからご応募下さいませ。アルバイトに関しては3ヶ月位の短期でも条件次第では可能ですので、お気軽にご相談下さいませ。


青木健児

データ分析部の青木です。専門は多変量解析です。階層ベイズモデルにも興味あります。普段の業務ではSQLとPython(pandas)をよく使います。