2012/09/04

文系のための「主成分分析の仕組み」(1)

主成分分析は、我々、文系研究者が最も良く利用する手法の一つである。
その一方で、誤用や誤解釈が多く、間違った使い方も多い。
そこで、本ブログでは、主成分分析の仕組みについて整理を行う。

そもそも、主成分分析とは何なのか?データを扱いながら考えてみる。
今回のデータも Google Spreadsheet に置いてあるので、
最初にRCurl ライブラリを読み込む。

library(RCurl)

データの詳細は、散布図行列の話にあるので、そちらを参照する。
ライブラリを無事に読み込めたら、次は、以下のコマンドを実行する。

# Google Spreadsheet からデータをダウンロードする
data <- getURL("https://docs.google.com/spreadsheet/pub?key=0AtOtIs5BjRVhdHo0c0pjb29OTU9aV3BmQUFJUWJQa0E&single=true&gid=0&range=A1%3AK48&output=csv")

# ダウンロードしたデータを読み込む
X <- as.matrix(read.table(textConnection(data), header=TRUE, row.names=1, sep=","))

データを扱う前に、最初にするべきことは何だったか?
確か、相関係数行列を見て、対応する散布図行列を参照するのだった。
# 散布図行列を作成するためにpsych ライブラリを読み込む。
library(psych)

# 相関係数行列を出す
round(cor(X), 3)

#  散布図行列を出す。
pairs.panels(X, scale=TRUE, smooth=FALSE, ellipses=FALSE, density=FALSE)

なお、相関係数は以下の通り。小数点以下が長いので、round()関数で丸め込み。

> round(cor(X), 3)
        HM    ELD     HC    CHR     SC     LI     SP    ENV    DIS     IC
HM   1.000  0.023  0.103 -0.106  0.201  0.124 -0.205  0.140 -0.134 -0.143
ELD  0.023  1.000  0.350  0.470 -0.046  0.418  0.158  0.341 -0.012  0.429
HC   0.103  0.350  1.000  0.310  0.082  0.275  0.077  0.162 -0.007  0.381
CHR -0.106  0.470  0.310  1.000 -0.001  0.564  0.495  0.430  0.055  0.182
SC   0.201 -0.046  0.082 -0.001  1.000 -0.082  0.001  0.027  0.072  0.062
LI   0.124  0.418  0.275  0.564 -0.082  1.000  0.323  0.754  0.097  0.007
SP  -0.205  0.158  0.077  0.495  0.001  0.323  1.000  0.341  0.121  0.170
ENV  0.140  0.341  0.162  0.430  0.027  0.754  0.341  1.000  0.044 -0.020
DIS -0.134 -0.012 -0.007  0.055  0.072  0.097  0.121  0.044  1.000  0.271
IC  -0.143  0.429  0.381  0.182  0.062  0.007  0.170 -0.020  0.271  1.000


さて、この相関係数行列散布図行列を見ることには慣れてきているはず。
そのように信じたいし、それを前提に話をするのであるが…。

ふむ。これでも、データの全体像を掴むことはできるのであるが、
あっちを見ながら、こっちを見て…」というのは、少々、面倒である。
せめて、散布図行列対角成分に入っているヒストグラムを見るだけなら…。

主成分分析というのは、要するにそういうことである。

元のデータを「対象のみの情報」と「変数のみの情報」に分け、
さらに、本質を変えることなく、対象のみの情報を観察する

解ったような、解らないような話であるが…とにかくやってみよう。

まずは、多次元データの分解の話からはじめる。
うん?「分解」って…どこかで聞いたような?

そうそう。その話は、入門編特異値分解の話で済ませてある。
実は、この特異値分解こそが、主成分分析を理解する上での「」となる。

今回は、珍しく、図を使って整理してみる。


普段、我々が扱っている多次元データというのは、
n個の対象」と「p個の変数」から成る「n×p行列」の「」である。
そして、個々の値のバラツキ(偏差)が集まった行列を「」と言った。

実は、この「」というのは、理想的というか観念的な世界では、

対象のみの本質であるn×n行列」と、
変数のみの本質であるp×p行列」として存在している。

この二つの行列が上手く組み合わさったものがである。

観念論を受け入れることができる文系研究者には解りやすいだろう。

さて、結局のところ、特異値分解における二つの特異値ベクトルは、
一方が、「対象の本質」を表した行列の「回転装置」で、
もう一方が、「変数の本質」を表した行列の「回転装置」である。

ここで出てきた二つの行列は「歯車のようなもの」であるが、
二つの歯車の「ギア比」は、当然のことなが、異なっていてかみ合わない。
特異値は、この二つの歯車が上手く噛みあうような歯車と考えることができる。

主成分分析を行うと、「主成分得点」と「主成分負荷」という言葉が登場するが、要するに、

変換が掛かった行列Uギア比を調整するための「Σ」を掛けたものが主成分得点であり、
変換が掛かった行列Vギア比を調整するための「Σ」を掛けたものが主成分負荷である。

なるほど。実は、特異値分解主成分分析というのは、密接な関係があって、
特異値分解を通して主成分分析の仕組みを考えると非常に解りやすい。

ところで、この特異値「Σ」というのは、面白い性質を持っている。
実は、特異値は、大きい方から並んで出てくる性質があって、
大きい値ほどギア比の調整に大きな影響を持っている。

どういうことか?先ほど、ギア比で説明したが、
二つのギア比を「一回の調整」だけで完全に合わせることは非常に困難で、
何回かの調整を追加で行なって徐々に調整していく必要がある。

精度を上げるということは、何回も微調整を続けることで、
つまりは、この調整作業は、ベクトルとして出てくる特異値を。
前から順に、何個目まで使うか?ということに他ならない。

見方を変えると、最後の方の特異値は「微細な」補正程度のものなので、
前の方のいくつかの特異値を使うだけで、データの構造はある程度復元できてしまう

主成分分析では、「第何主成分まで採るか?」という問題があるが、
要するに、「復元精度をどの程度で見切りをつけるか?」という話なのである。

試しに、読み込んだデータを用いてやってみよう。
まずは、偏差行列を計算するところから。

# ステップ1:データの総数(n)を求める。
n <- nrow(X)
# ステップ2:卸売数量の平均を求める。転置が必要。
m <- t(colMeans(X))
# ステップ3:行列の引き算ができるように、平均値を総数個複製する。
M <- matrix(rep(m, each=n), nrow=n)
# ステップ4:元の行列から、偏差行列を作成する。
A <- as.matrix(X-M)

次に、偏差行列が入った行列を特異値分解し、A.svdという変数に格納する。
特異値分解は、確か、svd()関数で出来た。
さらに、特異値分解した結果を用いて、主成分得点主成分負荷量を計算する。

# 特異値分解を行う
A.svd <- svd(A)

# 主成分スコアと主成分負荷量を計算する
A.svd.score <- A.svd$u %*% diag(A.svd$d)
A.svd.loadings <- A.svd$v %*% diag(A.svd$d)

# 結果を出力する
head(A.svd.score)
head(A.svd.loadings)

演算結果は以下の通り。全てを出力すると大変なので、一部のみ。

> head(A.svd.score)
           [,1]      [,2]      [,3]       [,4]       [,5]      [,6]      [,7]
[1,]   3.651314 -7.957318  5.140459   1.052492  2.3098159 -4.752406  1.800883
[2,]  19.786392  5.057985 12.887751 -11.766562  0.6327105  3.429392  5.962028
[3,]   4.191283 18.645826 -4.698321   7.407312  3.6354721  2.194497  2.183438
[4,] -17.587890 -1.574478  2.599470  -2.081724 -9.3177502 -1.848632 -3.385567
[5,]  -4.172353 16.769956 10.395899   3.748746  6.7418704  5.452591  6.658994
[6,]  21.314091 14.733086  1.864985   9.509779  1.0878435 -2.527682  2.741948
          [,8]       [,9]       [,10]
[1,] 0.2884762  0.7450070  1.88445215
[2,] 0.7499275 -1.0444792 -0.18379527
[3,] 3.0996914 -2.9885232  1.41489124
[4,] 6.7199681 -5.7687203 -1.24307489
[5,] 3.5853302  0.8062032  0.03691304
[6,] 0.8179828  0.5430898  1.97413603

> head(A.svd.loadings)
           [,1]       [,2]       [,3]        [,4]        [,5]       [,6]
[1,]   4.091736  -9.947145   6.385787 -10.9269664 -26.1646501   9.507016
[2,] -21.716695 -15.885216 -13.466549   2.7171999 -12.7217615 -26.129684
[3,] -37.546146 -45.813350  32.850625  12.0989336   5.1090391   2.096065
[4,]  -6.881634 -12.306043  -8.587840   0.3087546   6.7732727  -7.116862
[5,]  -3.811058  -3.656824  15.081127 -45.3911883   6.3497666  -6.610114
[6,]  -2.508814 -16.633645 -11.695187  -0.2798455  -0.4216687   1.218834
              [,7]       [,8]       [,9]      [,10]
[1,] -14.361831207  1.1456268 -0.1570407 -0.7763315
[2,]   0.134726725 -2.0857931 -4.1224558 -0.3527740
[3,]   0.968825495 -0.6862824 -1.1022974 -0.3394291
[4,]  -7.044272819  9.5553820 10.2820549 -4.4587931
[5,]   2.353135513 -0.1068618  0.2054425  0.7201806
[6,]  -0.009821734  7.0229551  1.2055759 10.5304434

以上のようにして、主成分得点主成分負荷量を無事に出すことができた。

ところで、冒頭で「分散の部分だけ見れれば…」という話をした。
これはどういう意味か?試しに、主成分得点分散を見てみる。

# 小数点が長すぎるので、小数点第三位までで丸め込み。
round(var(A.svd.score),3)

実行結果は以下の通り。「非常に面白いこと」になっている。

> round(var(A.svd.score),3)
         [,1]    [,2]   [,3]   [,4]   [,5]   [,6]   [,7]  [,8]  [,9] [,10]
 [1,] 312.518   0.000  0.000  0.000  0.000  0.000  0.000 0.000 0.000 0.000
 [2,]   0.000 100.764  0.000  0.000  0.000  0.000  0.000 0.000 0.000 0.000
 [3,]   0.000   0.000 65.603  0.000  0.000  0.000  0.000 0.000 0.000 0.000
 [4,]   0.000   0.000  0.000 53.541  0.000  0.000  0.000 0.000 0.000 0.000
 [5,]   0.000   0.000  0.000  0.000 28.354  0.000  0.000 0.000 0.000 0.000
 [6,]   0.000   0.000  0.000  0.000  0.000 22.807  0.000 0.000 0.000 0.000
 [7,]   0.000   0.000  0.000  0.000  0.000  0.000 15.126 0.000 0.000 0.000
 [8,]   0.000   0.000  0.000  0.000  0.000  0.000  0.000 7.267 0.000 0.000
 [9,]   0.000   0.000  0.000  0.000  0.000  0.000  0.000 0.000 6.065 0.000
[10,]   0.000   0.000  0.000  0.000  0.000  0.000  0.000 0.000 0.000 3.257

なんと、主成分得点というのは、変数の関係から解き放たれているので、
変数同士の関係が無い状態、つまり、「無相関」という状態になる。
したがって、共分散が全て「0」となり、分散のみが対角成分に入るのである。

つまり、今までは、変数間の関係に気を配りながら散布図行列を観察していたが、
対角線上に並ぶヒストグラムのみに気を配れば良い状態になる。

ここで、もう一度、ギア比の話を思い出してみる。
そうそう、「調整を何回するのか?」という話。

実は、この問題もあっさり解決できる。

どういうことか?今度は、主成分得点の分散の比率を見てみる。

# 主成分得点の分散を変数に格納する
A.svd.score.var <- diag(var(A.svd.score))

# 主成分得点の分散の合計を変数に格納する。
A.svd.score.var.sum <- sum(A.svd.score.var)

# 主成分ごとの分散の割合を計算して表示する
for(i in seq(1:length(A.svd.score.var))){
    print((A.svd.score.var[i]/A.svd.score.var.sum)*100)
}

以下は、その実行結果。

> # 主成分ごとの分散の割合を計算して表示する
> for(i in seq(1:length(A.svd.score.var))){
+     print((A.svd.score.var[i]/A.svd.score.var.sum)*100)
+ }
[1] 50.79092
[1] 16.37642
[1] 10.66193
[1] 8.701498
[1] 4.608147
[1] 3.706659
[1] 2.458381
[1] 1.181115
[1] 0.9856477
[1] 0.529284

実は、この分散の割合というのは、元の対象間の関係の説明の割合をしめしている。
ギア比」の説明の通り、特異値最初から何個まで採るかが、精度を決めると書いた。

したがって、この結果から上から順に、
  • 第一主成分の説明力:50.7%
  • 第二主成分の説明力:16.4%
  • 第三主成分の説明力:10.7%
  • 第四主成分の説明力:8.7%
  • 第五主成分の説明力:4.6%
  • 第六主成分の説明力:3.7%
  • 第七主成分の説明力:2.5%
  • 第八主成分の説明力:1.2%
  • 第九主成分の説明力:1.0%
  • 第十主成分の説明力:0.5%
となる。さらに、上記の式を少し変えて、累積比率を出してみると、

# 累積比率の変数を初期化する
vsum <- 0

# 主成分ごとの分散の累積割合を計算して表示する
for(i in seq(1:length(A.svd.score.var))){
    vsum = vsum + (A.svd.score.var[i]/A.svd.score.var.sum)*100
    print(vsum)
}

となり、その結果は、以下の通り。

> # 主成分ごとの分散の累積割合を計算して表示する
> for(i in seq(1:length(A.svd.score.var))){
+     vsum = vsum + (A.svd.score.var[i]/A.svd.score.var.sum)*100
+     print(vsum)
+ }
[1] 50.79092
[1] 67.16734
[1] 77.82927
[1] 86.53077
[1] 91.13891
[1] 94.84557
[1] 97.30395
[1] 98.48507
[1] 99.47072
[1] 100

つまり、第四主成分までで、元のデータの86%が説明できる
考えようによっては、第五主成分以下を切り捨てても良いかもしれない

この状況を上手く、プロットしたものが主成分プロットなのだけれども、
その話は、今回はしない。まだまだ、基本的な話は終わっていない。

さて、実を言うと、この主成分得点の分散というのは、
分散共分散行列の特異値分解の特異値に相当する。つまり、

svd(cov(X))$d

ということで、実際に実行すると以下のようになる。

> svd(cov(X))$d
 [1] 312.517713 100.764478  65.603086  53.540516  28.354034  22.807157
 [7]  15.126472   7.267429   6.064713   3.256697

なるほど、つまり、分散共分散行列特異値分解を行うと、
その特異値は、主成分得点の分散に一致し、したがって、
元の行列から分散だけ見れば良い状況が得られるである。

2 件のコメント:

  1. 「# 主成分得点の分散の合計をへんすに格納する。」
    変換ミスの様ですが・・・?

    返信削除
    返信
    1. 変換ミスですね...修正しておきました。
      コメントありがとうございます。

      削除