文系のための「主成分の選択」

主成分分析の仕組みも理解できた。解釈の仕方も理解できた。
これで、基本的に主成分分析を使うことは出来る。

しかしながら、当初に問題として挙げておいて、まだ整理できていない問題がある。
すなわち、主成分を何番目まで採るべきか？という問題である。
これには、いくつかの方法が存在し、ここでは、4つの方法について整理する。

とりあえずは、これまでのデータを準備する。

library(RCurl)

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

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

この結果を主成分分析にかけて、さらに、データの要約を確認する。

X.pca <- prcomp(X)
summary(X.pca)

二行に分かれてしまっているが、とにかく、以下が実行結果。

> X.pca <- prcomp(X)
> summary(X.pca)
Importance of components:
                           PC1     PC2    PC3     PC4     PC5     PC6     PC7
Standard deviation     17.6782 10.0382 8.0996 7.31714 5.32485 4.77568 3.88928
Proportion of Variance  0.5079  0.1638 0.1066 0.08701 0.04608 0.03707 0.02458
Cumulative Proportion   0.5079  0.6717 0.7783 0.86531 0.91139 0.94846 0.97304
                           PC8     PC9    PC10
Standard deviation     2.69582 2.46266 1.80463
Proportion of Variance 0.01181 0.00986 0.00529
Cumulative Proportion  0.98485 0.99471 1.00000

Standard deviation は分散、Propotion of Variance は分散の割合、
そして、Cummulative Proportion は累積された分散の割合を示す。
特に、最後の累積された分散の割合のことを累積寄与率と呼ぶ。

さて、分散の割合がなぜ、重要であったかを思い出してみる。
たしか、主成分分析のしくみの話でも述べたこと。
つまり、元のデータの説明力を示しているのであった。

主成分の採用基準の1つめの方法というのは、この累積寄与率を見る方法である。
よく見かける方法は、「大体、7割程度の近似で良いのではないか？」というもの。

したがって、この場合には、第3主成分までを採用すると、
累積寄与率が、約78%となり、第3主成分まで判断すれば良いことになる。

しかしながら、この基準の場合、困ったことがある。
非常に、累積寄与率が7割付近で、微妙に推移する場合、
その微妙な差で、採用と不採用を決めることができるか、という問題である。

では、二つ目の方法。

今度は、スクリープロットというものを用いる。

# スクリープロットを描画する。
screeplot(X.pca, main="Screeplot for Variance")

スクリーとは、「崖」を意味し、急激な変化が起きている部分までを採用するという方法。
Rでは、screeplot()関数を用いる。それにしても、何とも、愛想の無いプロットである。

この場合は、第1と第2主成分の間、第4と第5主成分の間に変化が見える。
したがって、第1主成分のみ、あるいは第4主成分までを採用するということになる。

ふむ。もっと、スッキリとした方法は存在しないのであろうか？
ということで、三つ目の方法。

今度は、もう少し、明確な基準。カイゼル基準と呼ぶ。

この方法は、相関係数行列に対応する主成分分析を行い、
分散が「1」以上となるものを採用するというもの。

すでに述べたように、相関係数行列に対応する行列による主成分分析の場合、
主成分得点の分散の平均は「1」となった。

そもそも、標準化した行列による主成分分析であるので、
元の状態では、各変数における値の平均は「0」、分散は「1」である。そして、
主成分分析を行うと、変数間の関係から開放され、分散の合計が再配分される。

ここで視点を変えてみると、分散が「1」よりも小さいということは、
元の状態よりも、少ない状況しか反映できていないことを示す。
折角なので、先程のスクリープロットと合わせて表示してみる。

# とりあえず、標準化したデータによる主成分分析を行う
X.pca.scale <- prcomp(X, scale=TRUE)

# 一応、分析結果の要約も表示
summary(X.pca.scale)

# スクリープロットを描く
screeplot(X.pca.scale, main="Screeplot for Variance (with scaled data)")

# 分散が「1」となる基準線を弾く。
abline(h=1, lty=2)

要約の結果は以下の通り。

> # 一応、分析結果の要約も表示
> summary(X.pca.scale)
Importance of components:
                          PC1    PC2    PC3    PC4     PC5     PC6     PC7
Standard deviation     1.7355 1.2277 1.1613 1.0490 0.93393 0.79329 0.75419
Proportion of Variance 0.3012 0.1507 0.1348 0.1100 0.08722 0.06293 0.05688
Cumulative Proportion  0.3012 0.4519 0.5868 0.6968 0.78403 0.84696 0.90384
                           PC8     PC9    PC10
Standard deviation     0.68435 0.54333 0.44506
Proportion of Variance 0.04683 0.02952 0.01981
Cumulative Proportion  0.95067 0.98019 1.00000

第5主成分以下は、分散が「1」よりも低いので、第4主成分までが採用となる。
なるほど。この方法は、かなり良い。先程までの方法と比較しても客観的に見える。

では、この方法が最良と言えるかと言うと、ここにも問題がある。
主成分得点の分散は、後ろに行けば行くほど、値が小さくなるのは当たり前。
なんとか、この部分の調整をしたい。

ということで、四番目の方法、「平行分析」。「原理上」は、これがベスト。

この方法は、元のデータサイズと同じサイズ（対象と変数の数が同じ）の、
人工的な架空の値が入った行列を大量に準備し、そのデータから基準を作る方法。

一般的な方法は、人工的な架空の値には、標準正規乱数というのを用いる。
ここでは、標準正規乱数に関する詳しい話は省略するが、
要するに、誤差が平均「0」分散「1」となるような乱数を発生させる。

このようにして、元のデータと同じサイズの行列の成分を乱数で埋め尽くし、
そのようにして作られた行列で、主成分得点の分散を求める。

ここで、各々の成分は、乱数なので、主成分分析を1回だけ実行しただけだと、
誤差が影響してしまうので、それを1000回、10000回というように何回も実行する。

そうして、繰り返し計算で得られた主成分得点の分散の平均値、中央値、あるいはIQRと、
本来の分析対象であるデータの主成分得点の分散を比較し、上側に来るものを採用する。
要するに、カイゼル基準の改良版といったところであろうか。

なんとも、複雑そうに見えるが、Rでは、この方法は簡単に実行できる。
とは言え、標準では入っていない。psych() ライブラリを読み込む。

library(psych)
fa.parallel(X)

そして、以下が実行結果。

> fa.parallel(X)
Loading required package: MASS
Parallel analysis suggests that the number of factors =  3  and the number of components =  1

psych()関数のfa.parallel()関数を用いると、主成分分析の主成分の採用基準と、
同様に、因子分析の因子数の採用基準が表示される。

この図において、✕印の折れ線が主成分分析の場合を示し、
△印の折れ線が因子分析の場合を示している。

赤色の点線は、標準正規乱数による基準の線で、
破線は、別の方法（元のデータの入れ替えの方法）による基準線を示している。

出力結果を見てみると、以下のように書いてある。

「平行分析は因子=3、および主成分=1を提案する」

平行分析は、因子分析にも対応しているため、両方の結果を提案してくれるのだが、
今回は、主成分分析なので主成分について観察する。

さて、今回の場合、結果として主成分は「1」であるとのお告げがあったのだが、
これは平均値を基準にしているので、多少の余裕（IQR程度）を見ても構わない。
そこで、平行分析のプロットをもう一度観察してみる。

どうやら、第2主成分はごく僅かに基準線よりも下側に有り、
第3主成分は基準線よりも僅かに上側に来ている。
したがって、微妙ではあるが、第3主成分までは採れそうである。第4主成分以下は不可。

さて、最後に重要なこととして、主成分の採用基準というのは、
あくまで、元のデータの説明力に意味があるという前提がある。
一方で、重要な情報が元の情報の多さと一致しているとは限らない。

つまり、主成分の分散の大きさに目が行き過ぎるのも問題で、
実は、主成分の後ろの方に重要な意味が潜んでいる場合もある。

したがって、第1主成分と第2主成分で上手く説明できない場合には、
寄与率が低くとも、様々な組み合わせを確認することも重要である。

文系のための「主成分分析の可視化」（２）

主成分分析の解釈の仕方が理解できたところで、
いよいよ、一般的な主成分分析の可視化方法である「バイプロット」を実行してみる。

今回の話は、ソフトウェアが既定で設定している値を鵜呑みにしてはいけない例。
数理的な背景を理解していれば、簡単に理解できることである。

「文系だから深く知る必要は無い。使えれば良い」などと言う人がいる。
学生であっても、時には研究者もそのようなこと言う。最低限の知識は必要。
何度も主張しているが、私は、数学が苦手な文系研究者である。

さて、「最低限の知識」との線引きは確かに困難であるが、
少なくとも、間違った解釈をしないための知識は必要であろう。
特異値分解は、必要最低限の知識に加えても良いと思う。

特異値分解を概念的に理解しておけば、逆行列や主成分分析の理解も楽である。

一般的な教科書では、主成分分析を固有ベクトルと固有値から説明するが、
文系人間には、理解しにくく、何より、実際の計算では特異値分解を使っているので、
やはり、特異値分解から説明した方が良いであろう。

話が逸れたが、さっそく、バイプロットを使って、主成分分析の結果を可視化する。
本来ならば、分析前に散布図行列が必要であるが今回は省略。
データの詳細は、散布図行列の話にあるので、そちらを参照する。

library(RCurl)

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

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

まずは、prcomp()関数 を用いて、バイプロットを作成する。
今回は、標準化を掛けない方法と、掛けたものを比較する。
それぞれ、分散共分散行列と相関係数行列に対応するのだった。

# 並べて描くためのオマジナイ
par(mfrow=c(1,2))

# 標準化を行わない場合
biplot(prcomp(X), main="PCA without scaling")

# 標準化を行う場合
biplot(prcomp(X, scale=TRUE), main="PCA with scaling")

このように、二つの図を重ね合わせることで、より解釈がしやすくなる。

なお、バイプロットにおける二つの図の軸スケールは関係せず、
変数間の主成分負荷量を表す矢印線の向きおよび長さの違いと、
その方向に布置されている各対象の布置されている状況のみが重要である。

さて、Rにおけるバイプロットは、biplot()関数を使う。
基本的には、この方法でプロットすることが可能であるが、
既定の状態では、いくつか問題がある。

そもそも、バイプロットによって重ね合わされた二つの図、つまり、
主成分得点と主成分負荷量は根本的に異なる空間に存在している。
それを、半ば「強引に」重ね合わしている。

それ故に、重ねあわせ方についてはいくつかの主張が存在する。
その主張とは、すなわち、特異値の掛ける比率であり、
この比率は、二つの図を重ねあわせる際の回転に影響する。

タイプ	得点	負荷
Type 1	U✕D	V✕D
Type 2	U	V✕D
Type 3	U✕√D	V✕√D
Type 4	U✕D	V
Type 5	√n✕U	V÷√n

ここでは、便宜上、5つの形式に分けた。
biplot()関数の既定では、Type 2 が採用されており、
PCA()関数では、Type 1 が採用されている。

さて、ここで、Type 5 を除く、Type 1 〜 Type 4 の主張は、
特異値分解と大きく関わりがある。



この式は、偏差行列に対する特異値分解である。

ここで、Type 1は、特異値Σ分だけ余分であるが同じ空間に投影できる。
一方、Type 2 〜 Type 4 までは、特異値分解に合わせた布置となり、
Type 2 と Type 4 を採用した場合、特異値を余分に掛けている分だけ角度は合わない。

実際に角度を合わせるのであれば、Type 1 か Type 3 を採用するべきである。
この二つの方法は、変数と対象の配置関係は同じであるが、軸スケールは異なる。

さて、このように状況を整理した上で、再び、prcomp()関数について整理してみる。
実は、R の biplot()関数では、Vに乗じるDの割合を指定することができ、
それは、「0〜1」の連続的な数値で設定することができる。

つまり、特異値分解を以下のように書き換えた状況で考えている。



この調整分を行うパラメータは、 scale であり、既定では「1」となっている。
つまり、Uに対してはDを掛けずに、VのみにDを掛けている状態（Type 2）となる。
そして、この値を「0.5」とすると、Type 3 の方式となる。

# 並べて描くためのオマジナイ
par(mfrow=c(1,2))

# 標準化を行わない場合
biplot(prcomp(X), main="PCA without scaling", scale=0.5)

# 標準化を行う場合
biplot(prcomp(X, scale=TRUE), main="PCA with scaling", scale=0.5)


先程の図と見比べると、全体的に右方向に傾いていることが解る。
これは、特に、矢印の長いもので確認すると解かりやすい。

バイプロットは、変数の状況と対象の状況を比較して解釈する上で重要であるが、
微妙な角度の違いによって、間違った解釈を導きかねない。
したがって、これらの数理的背景を含めて理解することが重要である。

文系のための「主成分分析の可視化」（１）

前回までは、主成分分析の仕組みについて整理してきた。
主成分分析というのは、分散共分散行列と相関係数行列と深い関係があって、
これらの見通しをより良くするために、ある種の変換を掛けるであった。

重要なことは、変数が多いと変数間の関係が複雑になりすぎるので、
変数間の関係を本質を変えずに取り除き、対象間の関係のみに焦点を絞ること。

ところで、主成分分析を2〜3変数のデータに対して用いている例を見かけるが、
実は、少ない変数の場合は、分散共分散行列と相関係数行列を見ても混乱しない。

したがって、散布図行列を観察するだけで十分に解釈できる場合が多い。
状況に合わせて、適切な方法を選択することが重要である。
重回帰分析や偏相関係数を観察するという方法もある。

さて、少し話が逸れたが、本日は、主成分分析の解釈について整理する。
正確には、主成分得点と主成分負荷量を可視化し、
それらから適切な解釈を導く方法について整理する。

では、さっそく、前回までに用いたデータを用いて、実際の処理を始める。
今回もGoogle Spreadsheet にデータを置いてあるので、それを利用する。
なお、今回は、散布図行列を作成するステップは省略する。本当は必要。

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=","))

主成分分析の結果を可視化する場合には、
バイプロットと呼ばれる方法を用いるのが一般的であり、
prcomp()関数の結果も、その方法で可視化する。

しかしながら、最初からバイプロットで可視化すると、
初学者の多くは、間違った解釈を行う危険性がある。
そこで、本日は、特異値分解から直接的に分析し、その結果の可視化を行う。

ところで、このチュートリアルで用いているデータは、
分散共分散行列と相関係数行列のうち、どちらが適切であったか？

確か、単位が揃っていて、その単位における値の大小にも意味があったので、
分散共分散行列の方が適切であった。前回の話では、そのように解説した。
っということは、つまりは、偏差行列を準備する必要があった。

# 行列の引き算ができるように、平均値を総数個複製する。
M <- matrix(rep(m <- t(colMeans(X)), each=nrow(X)), nrow=nrow(X))

# 元の行列から、偏差行列を作成する。
A <- as.matrix(X-M)

次に、作成された偏差行列に対して特異値分解を行い、
主成分得点と主成分負荷量を計算する。

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

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

各主成分の説明力についての情報も欲しい。説明力は何に対応していたか？
確か、主成分得点の行列の分散の割合に相当していた。

# 割合を100分率で表し、さらに、小数点第二位で切り詰める。
A.svd.per <- round(diag((var(A.svd.score)/sum(diag(var(A.svd.score)))))*100,2)

割合なので、小数点以下をダラダラと書いても仕方がない。
小数点第二位〜第三位当たりで切る。この場合は、第二位までで十分。

では、さっそく、「第1主成分」における「主成分負荷量」を可視化してみる。
今は、量の大きさを可視化したいので、棒グラフを用いる。詳しくは、棒グラフの話。

# 第一主成分というのは、主成分の第1次元ということ
dim <- 1

# 各変数の分散の割合の大きさを表すだけなので、とりあえず棒グラフで表現
barplot(as.vector(A.svd.loadings[,dim]), names=colnames(X))

# タイトルは必ず必要。タイトル文字を作成する。少々面倒。
pca.loadings.title <- paste("Factors of principal component in dimension", dim, sep="")
pca.loadings.title <- paste(pca.loadings.title , A.svd.per[1],sep="\n")
pca.loadings.title <- paste(pca.loadings.title , "%", sep="")

# タイトルをグラフに与える。
title(main=pca.loadings.title, ylab="variance")

まずは、主成分負荷量の見方から整理する。

この図では、第1主成分における主成分負荷量は、正と負の値を取っている。
実は、主成分分析の正と負は計算方法に依存しソフトや関数によって逆転する。

したがって、分散の大きさの差と、「0」を中心に二つの相反する性質を表す。
正と負は関係無いので、正負に捕らわれずに判断することが重要。

この図の場合、第1主成分では、国際協力に関する活動と医療や健康に関する活動、
の二つの変数に相当する主成分負荷量が両極にあって、国際協力に近い方から見ていくと、
障害者に関する活動、高齢者に関する活動・・・という順番で並んでいる。

医療や健康に関する活動「のみ」が反対側に向いていて、
しかも、出ている度合いがかなり小さいのが気になるが...とりあえず、次に進む。

こうした状況を観察してみると、第1主成分における二つの対極というのは、
一方がグローバルな活動を表し、もう一方がローカルな活動を表している「らしく」、
この状況は、元のデータの50.79％を説明していることが解る。

「らしく」というのは、この状況だけでは、まだ、なんとも言えない。焦ってはいけない。
主成分得点に関しても同様の図を作成し、
意味については、両方を比較しながら再度検討する。

# 第一主成分というのは、主成分の第1次元ということ
dim <- 1

# 各変数の分散の割合の大きさを表すだけなので、とりあえず棒グラフで表現
barplot(as.vector(A.svd.score[,dim]), names=rownames(X), las=2)

# タイトルは必ず必要。タイトル文字を作成する。少々面倒。
pca.loadings.title <- paste("Scores of principal component in dimension", dim, sep="")
pca.loadings.title <- paste(pca.loadings.title , A.svd.per[dim],sep="\n")
pca.loadings.title <- paste(pca.loadings.title , "%", sep="")

# タイトルをグラフに与える。
title(main=pca.loadings.title, ylab="variance")

第1主成分における主成分得点は、第1主成分における主成分負荷量に対応している。
つまり、この図においては、「0」を中心に下側はグローバルな活動の強さを表し、
一方、上側はローカルな活動の強さを表している。

なお、「0」付近というのは、どちらにも寄らないことを意味するので、
比較的に平均「的」な対象が集まっている。

そうそう、具体的な数値もあった方が良い。

> A.svd.score[(order(A.svd.score[,1])),][,1]
      Oita       Aichi       Miyagi     Ibaraki    Ishikawa    Kanagawa
-84.3659879 -27.2237009 -17.5878902 -17.4725893 -15.6203251 -15.3652856
       Gifu        Nara   Fukushima     Saitama       Ehime       Hyogo
-14.7809511 -12.8015618 -11.9976794 -11.9089262  -9.1302220  -8.3305146
   Kumamoto       Kochi      Toyama    Tokyo-to       Akita     Niigata
 -7.8976958  -6.4295937  -5.7279074  -4.9173176  -4.1723528  -3.8796902
   Shizuoka     Okinawa       Shiga    Osaka-fu     Fukuoka   Tokushima
 -2.8219903  -1.7026102  -1.2568486  -0.8619497   0.3182618   2.6431475
     Nagano       Gumma    Hokkaido      Iwate        Chiba     Okayama
  2.6834044   2.8809400   3.6513143   4.1912833   4.8209478   5.4525580
     Kagawa     Tochigi   Kagoshima   Hiroshima        Saga        Mie
  6.0702498   9.7733572   9.8101319  13.1450657  13.4666745  13.5595689
  Yamanashi     Tottori    Wakayama    Miyazaki    Kyoto-fu   Yamaguchi
 13.8047619  14.2913455  14.4812483  15.6424106  15.8493474  15.8798751
   Nagasaki     Aomori      Shimane    Yamagata       Fukui
 18.3468538  19.7863922  21.2399576  21.3140915  23.1504013

まず、グローバルな活動の強さを見ると、大分県が圧倒的に強く現れており、
次いで、愛知県、宮城県、茨城県、石川県、神奈川県・・・と続く。

一方、ローカルな活動の強さを見ると、福井県が最も値が大きく、
山形県、島根県、青森県、長崎県、山口県・・・と続く。

全体的な印象として、グローバルな活動は、比較的人口が集中している地域が多く、
ローカルな活動は、過疎化などが問題になっている地域が多いように思える。

ここで、気になるのは、グローバルな活動が大分県で高いこと。
このように「なんで？」という疑問は良く生じる。

そういった時に、元データを参照すると新しいことが見えてくる。
常に、元データと分析結果を交互に行き来することが重要。

> X
            HM  ELD   HC  CHR   SC   LI   SP  ENV  DIS    IC
Hokkaido  11.3 35.1 34.1 21.2 44.3 14.4 15.5 24.1  5.4  17.9
Aomori    15.2 21.0 22.4 13.3 55.2  7.9  8.4 19.9  8.3   8.0
Iwate      6.0 26.1 12.5 15.9 33.1 10.6 15.6 16.9 11.8  26.9
Miyagi    22.8 40.6 32.8 20.7 43.4 13.4 14.3 22.6 15.2  40.4
Akita      4.9 19.5 25.1 16.1 41.2  7.4  9.8 12.8 10.4  33.2
Yamagata   7.4 26.5 13.4 15.8 33.3  9.5 11.4 13.0  5.9   8.7
Fukushima 14.8 30.5 16.2 15.9 40.9 10.5 16.4 29.6 15.2  41.8
Ibaraki    7.9 29.6 38.6 21.4 63.1  9.5 29.6 26.9 10.2  38.5
Tochigi   16.2 36.2 24.3 20.4 51.2 12.1 14.5 24.9  4.7  14.5
Gumma     14.6 32.4 23.1 18.8 41.9  8.7 16.8 26.5  5.9  23.3
Saitama   10.3 36.3 21.8 18.3 34.5 15.0 24.7 28.9  8.5  38.3
Chiba     10.9 29.4 33.1 20.6 42.2 15.5 24.7 22.6  7.1  17.6
Tokyo-to  11.4 41.9 36.8 24.5 37.9 26.1 18.6 48.1  6.6  23.8
Kanagawa  11.9 31.4 26.0 24.6 44.7 16.0 23.7 42.2  5.4  40.7
Niigata   10.3 35.6 36.8 19.0 49.6  8.4 20.6 23.6  5.2  24.8
Toyama    16.8 25.1 28.3 15.6 37.8 11.5 16.5 19.2  9.0  32.7
Ishikawa   6.9 41.7 13.4 25.3 35.0 10.4 23.3 20.3  9.1  43.9
Fukui      7.5 29.3 12.1 17.4 44.2 10.2 15.1 19.6  4.4   5.8
Yamanashi 12.8 27.4 18.2 17.1 39.2  8.7 12.0 25.4  5.4  14.6
Nagano    12.4 24.9 31.8 16.6 29.9  9.3 14.1 18.5  4.9  22.9
Gifu      18.5 29.9 20.7 16.4 44.6  8.2 18.7 26.3  5.2  43.8
Shizuoka  10.6 32.7 18.5 17.8 50.8 11.1 12.8 21.1  5.4  30.9
Aichi      9.0 42.6 21.1 20.2 36.0 11.5 22.1 20.7  5.2  54.2
Mie       23.6 26.8 14.2 17.6 45.4 12.1 15.8 22.6  4.6  16.3
Shiga      9.5 28.0 15.8 22.3 40.9 10.1 13.3 20.7  4.6  31.1
Kyoto-fu  17.6 32.7 22.3 30.3 38.2 18.0 22.4 34.5  5.8   8.3
Osaka-fu  11.9 28.8 41.4 22.2 36.3 13.6 17.5 38.9  7.4  21.2
Hyogo      7.0 40.0 39.8 26.4 42.2 15.6 16.1 23.3  7.8  27.4
Nara      16.3 34.0 36.7 24.9 39.9 11.1 17.2 23.6  7.6  35.0
Wakayama   8.2 38.2 25.1 16.7 39.6 13.2 10.8 33.7  8.2   8.9
Tottori   12.6 31.6 18.3 15.9 41.9  9.2 12.7 23.7  5.1  13.2
Shimane   17.4 28.1 12.3 18.6 31.9  9.9 20.9 20.5  6.6   8.3
Okayama   17.9 39.6 26.9 23.0 32.7 11.9 12.7 27.5  4.6  18.2
Hiroshima 10.0 33.8 25.2 23.1 36.2 12.5 20.7 24.9  4.2  10.8
Yamaguchi 14.2 26.2 31.4 19.7 32.9 11.2 17.5 24.0  5.5   7.9
Tokushima 19.3 26.0 29.8 15.7 46.5  9.7 12.0 17.5  6.9  23.2
Kagawa    18.2 30.8 21.7 20.1 54.5 11.1 12.6 20.7  8.4  20.5
Ehime      9.0 44.3 26.2 18.7 37.5  8.0 13.4 22.5  7.3  33.1
Kochi     19.1 38.3 27.6 19.7 43.0 14.6 11.9 34.7  5.4  30.7
Fukuoka    9.6 29.7 30.5 27.9 44.0 14.4 18.6 24.0 13.3  22.8
Saga      18.9 31.5 27.8 20.7 52.2 10.2 20.0 17.4  5.2  10.2
Nagasaki  26.7 33.4 25.0 18.8 55.1 11.7 10.7 26.9  4.3   6.1
Kumamoto  21.5 38.1 56.9 17.3 30.4 12.5 13.2 20.0  4.2  23.5
Oita      13.4 44.5 58.9 23.8 47.7 12.8 14.5 21.2  6.7 102.6
Miyazaki  13.6 23.8 26.0 20.1 38.7 13.7 21.9 33.0  6.9   9.6
Kagoshima 21.2 26.6 22.2 14.1 40.3 11.1 10.7 18.7  4.7  18.4
Okinawa   21.7 39.9 26.4 22.4 46.6 22.4 20.7 46.0 10.6  24.3

なるほど、確かに、大分県では国際協力に関わるボランティア活動の割合が、
他の都道府県と比較して、群を抜いて高い。では、福井県はどうであるか？

うん？今度は、あまり、何かの値が突出しているようには見えない。

ついでに、データのサマリと箱ヒゲ図も確認。

> boxplot(X, main="Average Days for Participation in Volunteer Activities", ylab="Average days in a year")

> summary(X)
       HM             ELD              HC             CHR    
 Min.   : 4.90   Min.   :19.50   Min.   :12.10   Min.   :13.30
 1st Qu.: 9.80   1st Qu.:27.70   1st Qu.:20.90   1st Qu.:16.90
 Median :12.80   Median :31.50   Median :25.20   Median :19.70
 Mean   :13.80   Mean   :32.35   Mean   :26.59   Mean   :19.85
 3rd Qu.:17.75   3rd Qu.:37.20   3rd Qu.:31.60   3rd Qu.:22.25
 Max.   :26.70   Max.   :44.50   Max.   :58.90   Max.   :30.30
       SC              LI              SP             ENV    
 Min.   :29.90   Min.   : 7.40   Min.   : 8.40   Min.   :12.80
 1st Qu.:36.90   1st Qu.: 9.80   1st Qu.:12.75   1st Qu.:20.40
 Median :41.20   Median :11.20   Median :15.80   Median :23.60
 Mean   :41.89   Mean   :12.05   Mean   :16.53   Mean   :24.98
 3rd Qu.:45.05   3rd Qu.:13.50   3rd Qu.:20.30   3rd Qu.:26.90
 Max.   :63.10   Max.   :26.10   Max.   :29.60   Max.   :48.10
      DIS               IC      
 Min.   : 4.200   Min.   :  5.80
 1st Qu.: 5.200   1st Qu.: 13.85
 Median : 5.900   Median : 23.20
 Mean   : 7.028   Mean   : 25.08
 3rd Qu.: 8.250   3rd Qu.: 32.90
 Max.   :15.200   Max.   :102.60

さて、当初はローカルな活動の強さを表していたのかと思ったが、
福井県のデータは、全体的に各値が低く、特に、国際協力に関する活動は著しく低い。

この傾向は、福井県の次の山形県にも現れている。

したがって、ローカルというよりは「反」グローバル的とするか、
ボランティアに活動に対する総合的な取り組み度合いと見ても良かもしれない。

今度は、第2主成分に関しても同じことをしてみる。
今回は、二つの図を並べて同時に出力する。

#並べて描くためのオマジナイ。
par(mfrow=c(1,2))

# 今度は次元が2になる。
dim <- 2

# 一つ目の図には、主成分負荷量の図
# 各変数の分散の割合の大きさを表すだけなので、とりあえず棒グラフで表現
barplot(as.vector(A.svd.loadings[,dim]), names=colnames(X))

# タイトルは必ず必要。タイトル文字を作成する。少々面倒。
pca.loadings.title <- paste("Factors of principal component in dimension", dim, sep=" ")
pca.loadings.title <- paste(pca.loadings.title , A.svd.per[dim],sep="\n")
pca.loadings.title <- paste(pca.loadings.title , "%", sep="")

# タイトルをグラフに与える。
title(main=pca.loadings.title, ylab="variance")

# 次は、主成分得点の可視化
# 各変数の分散の割合の大きさを表すだけなので、とりあえず棒グラフで表現
barplot(as.vector(A.svd.score[,dim]), names=rownames(X), las=2)

# タイトルは必ず必要。タイトル文字を作成する。少々面倒。
pca.loadings.title <- paste("Scores of principal component in dimension", dim, sep=" ")
pca.loadings.title <- paste(pca.loadings.title , A.svd.per[dim],sep="\n")
pca.loadings.title <- paste(pca.loadings.title , "%", sep="")

# タイトルをグラフに与える。
title(main=pca.loadings.title, ylab="variance")

実際の値も表示する。これは、結果のみ。

> A.svd.score[(order(A.svd.score[,2])),][,2]
    Tokyo-to     Kumamoto     Osaka-fu      Okinawa     Kyoto-fu        Hyogo 
-25.48330834 -19.54068817 -17.89270947 -17.53265633 -11.60217588 -10.09716695 
    Wakayama     Hokkaido     Nagasaki     Miyazaki        Kochi      Niigata 
 -7.99181874  -7.95731799  -7.78153288  -7.29715992  -6.45035170  -6.32440501 
       Chiba     Kanagawa    Yamaguchi      Okayama      Ibaraki         Nara 
 -6.11855778  -6.08651180  -6.01139834  -5.64678869  -5.49745999  -4.43666871 
   Hiroshima      Fukuoka      Tochigi         Saga       Miyagi        Oita  
 -4.06168230  -3.71780644  -3.37521039  -2.60332055  -1.57447800   0.03078564 
       Gumma       Nagano    Tokushima        Ehime       Kagawa      Saitama 
  1.82577203   3.22769906   3.37501834   3.66652229   3.79919709   3.97604392 
     Tottori      Aomori     Yamanashi    Kagoshima       Toyama         Mie  
  4.70810331   5.05798505   5.39623139   6.31970748   6.62959774   6.96469789 
     Shimane        Fukui         Gifu     Shizuoka    Fukushima       Aichi  
  8.14448030   9.14451299   9.99399234  10.06750081  11.39446940  12.86385463 
       Shiga     Yamagata     Ishikawa        Akita       Iwate  
 13.29618683  14.73308634  15.04994751  16.76995612  18.64582591 

第2主成分は、元情報の16.38％の情報を持っている。
主成分負荷量を見ると、障害者に関する活動と国際協力に関する活動が対立している。
全体的に見てみると、恒常的な状況と特殊な状況に分かれているように見える。

一方、主成分得点を見てみると、下側に最も値が高いのが東京都で、
熊本県、大阪府、沖縄県、京都府・・・と続く。
逆に、上側に最も値が高いのは、岩手県で、秋田県、石川県、山形県・・・と続く。

それぞれの上位を見てみると、下側には西日本が多く含まれていて、
上側には東日本が多く含まれているようにも見える。
文化の違いだろうか？それとも、人口密度の問題だろうか？

研究として踏み込むのであれば、より詳細な検討を必要とする。
元データがどのようにして得られたのか、あるいは、外国人の居住人口の割合、
高齢者や障害者の人口比率、犯罪件数、環境汚染といった他の情報も必要とする。
少なくとも、今回のデータでは、これ以上の解釈はできない。

これは、主成分分析に限ったことでは無いが、データの背景を知ることは極めて重要。
また、ある分析から導かれるのが解釈ではなく、新しい疑問であることも多い。

分析結果を安易に鵜呑みにしたり、
強引に分析結果を結びつけることも厳禁！
こういった状況を目にすることは多い。

さて、話が逸れたが、主成分分析の結果を可視化するために、
棒グラフを用いて説明してきたが、これは、各主成分を1次元ごとに観察するため。
次は、二つの次元を同時に可視化する方法を行うが1次元毎に解釈するのは同じ。

では、第1主成分と第2主成分の主成分負荷量を同時に表現してみる。

# 第1主成分と第2主成分の次元を設定する
dim1 <- 1
dim2 <- 2

# 第1主成分と第2主成分の主成分負荷量をベクトルに入れる
x <- A.svd.loadings[,dim1]
y <- A.svd.loadings[,dim2]

# x軸の軸ラベル
label_x <- ""
label_x <- paste("Dim.", dim1)
label_x <- paste(label_x, "(")
label_x <- paste(label_x, A.svd.per[dim1], sep="")
label_x <- paste(label_x, "%)")

# y軸の軸ラベル
label_y <- ""
label_y <- paste("Dim.", dim2)
label_y <- paste(label_y, "(")
label_y <- paste(label_y, A.svd.per[dim2], sep="")
label_y <- paste(label_y, "%)")

# まずは、プロット領域を作成し、ラベルを配置する
plot(x, y, xlab=label_x, ylab=label_y, type="n", main="Loadings of principal component")
text(x, y, colnames(X), cex=0.8)

# 原点0からの方向を矢印で表現する
arrows(x0=0, y0=0, x1=x, y1=y, length=0.1, col="red")

# 最後に、原点軸を入れて整形する
abline(h=0, v=0, lty=2, col="lightgrey")

主成分負荷量の見方は、棒グラフで見た方法と同じ。要するに、一次元で観察する。
主成分負荷量というのは、あくまで、各次元における主成分の方向の強さなので、
幾何的なベクトルとして考えて矢印で表現する。

今度は、主成分得点の図も出してみる。

# 第1主成分と第2主成分の次元を設定する
dim1 <- 1
dim2 <- 2

# 第1主成分と第2主成分の主成分負荷量をベクトルに入れる
x <- A.svd.loadings[,dim1]
y <- A.svd.loadings[,dim2]

# x軸の軸ラベル
label_x <- ""
label_x <- paste("Dim.", dim1)
label_x <- paste(label_x, "(")
label_x <- paste(label_x, A.svd.per[dim1], sep="")
label_x <- paste(label_x, "%)")

# y軸の軸ラベル
label_y <- ""
label_y <- paste("Dim.", dim2)
label_y <- paste(label_y, "(")
label_y <- paste(label_y, A.svd.per[dim2], sep="")
label_y <- paste(label_y, "%)")

# まずは、プロット領域を作成し、ラベルを配置する
plot(x, y, xlab=label_x, ylab=label_y, type="n", main="Loadings of principal component")
text(x, y, colnames(X), cex=0.8)

# 原点0からの方向を矢印で表現する
arrows(x0=0, y0=0, x1=x, y1=y, length=0.1, col="red")

# 最後に、原点軸を入れて整形する
abline(h=0, v=0, lty=2, col="lightgrey")

# 第1主成分と第2主成分の主成分負荷量をベクトルに入れる
x <- A.svd.score[,dim1]
y <- A.svd.score[,dim2]

# x軸の軸ラベル
label_x <- ""
label_x <- paste("Dim.", dim1)
label_x <- paste(label_x, "(")
label_x <- paste(label_x, A.svd.per[dim1], sep="")
label_x <- paste(label_x, "%)")

# y軸の軸ラベル
label_y <- ""
label_y <- paste("Dim.", dim2)
label_y <- paste(label_y, "(")
label_y <- paste(label_y, A.svd.per[dim2], sep="")
label_y <- paste(label_y, "%)")
# まずは、プロット領域を作成し、ラベルを配置する
plot(x, y, xlab=label_x, ylab=label_y, type="n", main="Scores of principal component")
text(x+1, y+1, rownames(X), cex=0.8)

# 原点0からの方向を矢印で表現する
points(x, y, pch=16, col="blue")

# 最後に、原点軸を入れて整形する
abline(h=0, v=0, lty=2, col="lightgrey")

主成分得点は、散布図の「ように」可視化する。あくまで「ように」。
重要なことは、主成分負荷量と同様に、各次元ごとに1次元で解釈すること。

さて、ここでは第1主成分と第2主成分における主成分得点を同時に布置している。
したがって、比較的近い所に布置されている物同士は、
この二つの主成分のみにおいて似た性質を持っていると言える。

また、「0」付近に布置されるもの（今回は存在しない）は、
可視化された二つの軸において、平均「的」な性質を持っていることを意味する。

さて、最後に、相関係数行列に対応した場合でのプロットも作成してみる。
今回の場合、分散共産分散行列に対応した方法が良いのであるが練習だけ。

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

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

# 二つのプロットを並べて描くためのオマジナイ
par(mfrow=c(1,2))

# 第1主成分と第2主成分の次元を設定する
dim1 <- 1
dim2 <- 2

# x軸の軸ラベル
label_x <- ""
label_x <- paste("Dim.", dim1)
label_x <- paste(label_x, "(")
label_x <- paste(label_x, A.svd.per[dim1], sep="")
label_x <- paste(label_x, "%)")

# y軸の軸ラベル
label_y <- ""
label_y <- paste("Dim.", dim2)
label_y <- paste(label_y, "(")
label_y <- paste(label_y, A.svd.per[dim2], sep="")
label_y <- paste(label_y, "%)")

# 主成分負荷量のプロット
# 第1主成分と第2主成分の主成分負荷量をベクトルに入れる
x <- A.svd.loadings[,dim1]
y <- A.svd.loadings[,dim2]

# まずは、プロット領域を作成し、ラベルを配置する
plot(x, y, xlab=label_x, ylab=label_y, type="n", main="Loadings of principal component", xlim=c(-1,1),ylim=c(-1,1),asp=1)
text(x, y, colnames(X), cex=0.8)

# 中心に円を描く
x.circle <- seq(-1, 1, by = 0.01)
y.circle <- sqrt(1 - x.circle^2)
lines(x.circle, y = y.circle)
lines(x.circle, y = -y.circle)

# 原点0からの方向を矢印で表現する
arrows(x0=0, y0=0, x1=x, y1=y, length=0.1, col="red")

# 最後に、原点軸を入れて整形する
abline(h=0, v=0, lty=2, col="lightgrey")

# 主成分得点のプロット
# 第1主成分と第2主成分の主成分負荷量をベクトルに入れる
x <- A.svd.score[,dim1]
y <- A.svd.score[,dim2]

# まずは、プロット領域を作成し、ラベルを配置する
plot(x, y, xlab=label_x, ylab=label_y, type="n", main="Scores of principal component")
text(x+0.1, y+0.1, rownames(X), cex=0.8)

# 原点0からの方向を矢印で表現する
points(x, y, pch=16, col="blue")

# 最後に、原点軸を入れて整形する
abline(h=0, v=0, lty=2, col="lightgrey")

相関係数行列に対応する主成分分析を行う場合、
主成分負荷量は、sqrt(n-1) で基準化する必要がある。

そして、これを可視化する場合には、分散が1であることを示すために半径1の円を描く。
別の言い方をすると、主成分負荷量のプロットの中心に円が描かれている場合には、
その分析結果が相関係数行列に対応していることが解るのである。

さて、今回は第1主成分と第2主成分のみの可視化を行った。
これは、元のデータの67.17%を再現している。
これが十分であるかは、もう少し議論がある。

基本的には、この値が高いほど元のデータに近づくのであるが、
必ずしも、値が高ければ良いとも言えない。

というのは、第1主成分に特定の変数の値が強く反映され、
他の状況が上手く読み取れない場合が存在する。
そのような場合、第1主成分を除いて、第2主成分以下で分析することがある。

また、ある特殊な状況においては、第1主成分と第2主成分の主成分得点のプロットが、
「U」の字に布置されることがあって（馬蹄形問題）、そのような場合には、
第1主成分と第3主成分で観察し直す必要がある。この話はいずれ。

今回は、主成分負荷量と主成分得点を可視化する方法を実践し、
さらに、その結果を横に並べて表示した。
実は、この二つの図を重ね合わせる可視化方法があって、
その方法のことをバイプロットと呼ぶ。

次回は、そのバイプロットについて整理する。

文系のための「主成分分析の仕組み」（２）

前回は、特異値分解を通して主成分分析の仕組みについて整理した。
主成分分析は、「対象の本質」と「変数の本質」に分離し、
「対象のみに焦点」を当てて分析する方法だった。

そして、対象の本質のみを記述した部分を主成分得点と呼び、
変数の本質のみを記述した部分を主成分負荷量と呼んだ。

また、対象のみの状態にすることで、変数間の関係から開放され、
分散のみに注目するだけで良い状況になった。

前回の話では、特異値分解を用いて主成分分析を行ったが、
もちろん、Rには主成分分析の関数が存在している。しかも、色々。

そこで、今回は最初にRにおける主成分分析の関数を使ってみることにする。

今回もGoogle Spreadsheet にデータを置いてあるので、それを利用する。

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(cov(X), 3)

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

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

分散共分散行列と相関係数行列の出力結果は以下の通り。

> # 散布図行列を作成するためにpsych ライブラリを読み込む。
> library(psych)
>
> # 分散共分散行列を出す
> round(cov(X), 3)
         HM    ELD      HC    CHR     SC     LI     SP    ENV    DIS      IC
HM   27.370  0.750   5.470 -2.063  7.600  2.305 -5.029  5.586 -1.925 -12.463
ELD   0.750 38.669  22.176 10.859 -2.048  9.253  4.594 16.180 -0.199  44.576
HC    5.470 22.176 103.638 11.745  6.029  9.985  3.673 12.545 -0.205  64.775
CHR  -2.063 10.859  11.745 13.820 -0.021  7.465  8.619 12.185  0.559  11.328
SC    7.600 -2.048   6.029 -0.021 52.300 -2.108  0.040  1.468  1.434   7.519
LI    2.305  9.253   9.985  7.465 -2.108 12.677  5.382 20.473  0.950   0.394
SP   -5.029  4.594   3.673  8.619  0.040  5.382 21.932 12.175  1.550  13.332
ENV   5.586 16.180  12.545 12.185  1.468 20.473 12.175 58.090  0.927  -2.488
DIS  -1.925 -0.199  -0.205  0.559  1.434  0.950  1.550  0.927  7.540  12.452
IC  -12.463 44.576  64.775 11.328  7.519  0.394 13.332 -2.488 12.452 279.265
>
> # 相関係数行列を出す
> 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

前回は、svd()関数を使って主成分分析を行った。まずは、その復習。

偏差行列を出す処理は簡略化。まとめて実行。

# 行列の引き算ができるように、平均値を総数個複製する。

M <- matrix(rep(m <- t(colMeans(X)), each=nrow(X)), nrow=nrow(X))

# 元の行列から、偏差行列を作成する。

A <- as.matrix(X-M)

# 特異値分解を行う

A.svd <- svd(A)

# 主成分スコアと主成分負荷量を計算する

A.svd.score <- A.svd$u %*% diag(A.svd$d)

A.svd.loadings <- A.svd$v %*% diag(A.svd$d)

これで主成分分析は完了。主成分得点と主成分負荷量が出せた。

今回は、この結果とRの主成分分析の関数との比較するところから始める。

実は、Rの主成分分析の関数は一つではなく複数存在し、少しずつ異なる。

R の関数は、世界中の大多数の研究者、技術者、学生らによって開発されているので、

様々な手法が存在し得る。そのあたりが、商用製品と大きく異なるところ。

雑多という批判を受けることもあるが、中身がはっきりと見えるし、

問題があれば自分で修正もできる。自分で開発した関数を公開することもできる。

それゆえに、従来の商用製品よりも信頼性は高いと言われている。

さて、話は少し逸れたが、有名な主成分分析の関数は以下の通り。

• princomp()関数：stats パッケージの主成分分析の関数。そのまま使える。
• prcomp()関数：同じく、stats パッケージの主成分分析の関数。そのまま使える。
• PCA()関数：FactoMineR パッケージの主成分分析の関数。インストールが必要

では、これらの内、どれを使うべきか？難しい問題である。

princomp()関数は、色々と問題があって使うべきではない。

prcomp()関数は、理論通りの方法。ただし、バイプロットと呼ばれる図では注意が必要。

PCA()関数は、非常に多機能。ただし、色々と値が調整されているので注意が必要。

とりあえずは、最も無難そうなprcomp()関数を使うことにする。

# prcomp()関数で主成分分析を行い、変数に格納する。

X.pca <- prcomp(X)

# 分析結果の格納状態を確認する。

str(X.pca)

結果がどのように格納されているかというと以下の通り。

> str(X.pca)

List of 5

 $ sdev    : num [1:10] 17.68 10.04 8.1 7.32 5.32 ...

 $ rotation: num [1:10, 1:10] 0.0341 -0.1811 -0.3131 -0.0574 -0.0318 ...

  ..- attr(*, "dimnames")=List of 2

  .. ..$ : chr [1:10] "HM" "ELD" "HC" "CHR" ...

  .. ..$ : chr [1:10] "PC1" "PC2" "PC3" "PC4" ...

 $ center  : Named num [1:10] 13.8 32.3 26.6 19.8 41.9 ...

  ..- attr(*, "names")= chr [1:10] "HM" "ELD" "HC" "CHR" ...

 $ scale   : logi FALSE

 $ x       : num [1:47, 1:10] 3.65 19.79 4.19 -17.59 -4.17 ...

  ..- attr(*, "dimnames")=List of 2

  .. ..$ : chr [1:47] "Hokkaido" "Aomori-ken " "Iwate-ken " "Miyagi-ken" ...

  .. ..$ : chr [1:10] "PC1" "PC2" "PC3" "PC4" ...

 - attr(*, "class")= chr "prcomp"

では、特異値分解で行った方法と順番に比較する。

今回は、長くなるので出力だけ。出力が同じであることが確認できれば良い。

• $ sdev は、主成分得点の標準偏差。

> # prcomp()関数の$sdev の中身

> X.pca$sdev

 [1] 17.678171 10.038151  8.099573  7.317139  5.324851  4.775684  3.889277

 [8]  2.695817  2.462664  1.804632

> # 特異値分解の場合。エラーが出るがここでは無視。

> sd(A.svd.score)

 [1] 17.678171 10.038151  8.099573  7.317139  5.324851  4.775684  3.889277

 [8]  2.695817  2.462664  1.804632

• $ rotation は、特異値ベクトル「V」のこと。つまり、「変数の本質」の回転。

> # prcomp()関数の$rotation の中身、長いので一部のみ。

> head(round(X.pca$rotation, 3))

       PC1    PC2    PC3    PC4    PC5    PC6    PC7    PC8    PC9   PC10

HM   0.034 -0.146  0.116 -0.220 -0.724  0.294 -0.544  0.063 -0.009 -0.063

ELD -0.181 -0.233 -0.245  0.055 -0.352 -0.807  0.005 -0.114 -0.247 -0.029

HC  -0.313 -0.673  0.598  0.244  0.141  0.065  0.037 -0.038 -0.066 -0.028

CHR -0.057 -0.181 -0.156  0.006  0.188 -0.220 -0.267  0.523  0.616 -0.364

SC  -0.032 -0.054  0.275 -0.915  0.176 -0.204  0.089 -0.006  0.012  0.059

LI  -0.021 -0.244 -0.213 -0.006 -0.012  0.038  0.000  0.384  0.072  0.860

> # 特異値分解の場合。名前は入ってないが同じ。

> head(round(A.svd$v,3))

       [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  [,10]

[1,]  0.034 -0.146  0.116 -0.220 -0.724  0.294 -0.544  0.063 -0.009 -0.063

[2,] -0.181 -0.233 -0.245  0.055 -0.352 -0.807  0.005 -0.114 -0.247 -0.029

[3,] -0.313 -0.673  0.598  0.244  0.141  0.065  0.037 -0.038 -0.066 -0.028

[4,] -0.057 -0.181 -0.156  0.006  0.188 -0.220 -0.267  0.523  0.616 -0.364

[5,] -0.032 -0.054  0.275 -0.915  0.176 -0.204  0.089 -0.006  0.012  0.059

[6,] -0.021 -0.244 -0.213 -0.006 -0.012  0.038  0.000  0.384  0.072  0.860

• $ center というのは、元のデータの真ん中。列平均値。特異値分解とは無関係。

> # prcomp()関数の$center の中身。長いので一部のみ。
> round(X.pca$center ,3)
    HM    ELD     HC    CHR     SC     LI     SP    ENV    DIS     IC
13.804 32.349 26.585 19.849 41.885 12.053 16.532 24.983  7.028 25.081

> # 特異値分解とは無関係。元データの列平均。
> round(colMeans(X), 3)
    HM    ELD     HC    CHR     SC     LI     SP    ENV    DIS     IC
13.804 32.349 26.585 19.849 41.885 12.053 16.532 24.983  7.028 25.081

• $ scale なんだ、コレ!? 実は、これが今日の話のメインに関係。説明は後回し。

> # prcomp()関数の$scale の中身。とりあえず、中身の確認だけ。

> X.pca$scale

[1] FALSE

• $ x これは、主成分得点。

> # prcomp()関数の$x の中身。長いので一部のみ。
> head(round(X.pca$x ,3)[,1:5])
                 PC1    PC2    PC3     PC4    PC5
Hokkaido       3.651 -7.957  5.140   1.052  2.310
Aomori-ken    19.786  5.058 12.888 -11.767  0.633
Iwate-ken      4.191 18.646 -4.698   7.407  3.635
Miyagi-ken   -17.588 -1.574  2.599  -2.082 -9.318
Akita-ken     -4.172 16.770 10.396   3.749  6.742
Yamagata-ken  21.314 14.733  1.865   9.510  1.088

> # 特異値分解の場合。名前は入ってないが同じ。

> head(round(A.svd.score,3)[,1:5])
        [,1]   [,2]   [,3]    [,4]   [,5]
[1,]   3.651 -7.957  5.140   1.052  2.310
[2,]  19.786  5.058 12.888 -11.767  0.633
[3,]   4.191 18.646 -4.698   7.407  3.635
[4,] -17.588 -1.574  2.599  -2.082 -9.318
[5,]  -4.172 16.770 10.396   3.749  6.742
[6,]  21.314 14.733  1.865   9.510  1.088

以上のように、prcomp()関数は、特異値分解によって主成分分析をしていることが解る。
なお、主成分負荷量に関しては計算されておらず、
主成分プロットで負荷量をプロットする際に内部的に計算される。

さて、後回しにしていたが、唯一、なかったものが $ scale の部分だった。
これは一体、何だろうか？scale という言葉に見覚えは無いだろうか？
そうそう。確か、標準化のことを scaling と呼んだ。平均「0」で分散「1」。

実は、主成分分析は、分散共分散行列に対応させる方法と、
相関係数行列に対応させる方法の二種類の方法が存在する。

なぜか？よく考えてみると難しい話ではない。

要するに、各変数の単位が同じ場合と異なる場合があって、
単位が異なる場合には、全体的に大きな値の入った変数に影響されるから。

例えば、身長、体重、血圧、体温という変数があるデータがあったとする。
このとき、身長をミリ、体重をキロ、血圧をパスカル、体温を摂氏、
でそれぞれ計測していたとする。つまり、単位がバラバラの状況。

この状況というのは、実は、別の単位で置き換えることもできる。
身長をメートル、体重をグラム、血圧をミリ、体温を華氏、といった具合に。

この二つのデータは、単位が変わっただけで、本質的には何も変わらないが、
分散共分散行列の結果は、値に応じて変わってしまう。
一方、相関係数行列の場合は、平均と分散で調整しているため結果は同じ。

つまり、相関係数行列に対応するように主成分分析を行うことで、
単位の問題を解決することができる。

また、単位そのものが一緒であっても、変数間の値の差が明らかに大きい場合で、
かつ、そういった状況が望ましくないような場合には、
相関係数行列に対応した方法の方が適することがある。

さて、では、実際に相関係数行列に対応した方法を検討してみる。

# 偏差行列の代わりに、標準化したデータを用いる
Z <- scale(X)

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

# 主成分得点と主成分負荷量を変数に格納する

Z.svd.score <- Z.svd$u %*% diag(Z.svd$d)

Z.svd.loadings <- Z.svd$v %*% diag(Z.svd$d)

# 主成分得点の一部を表示する

head(round(Z.svd.loadings, 3)[,1:5])

この結果を、prcomp()関数と比較する。

# prcomp()関数のscaleパラメータをTRUEに変更する。

prcomp(X, scale=TRUE)

# 主成分得点の一部を表示する。

head(round(X.pca$x ,3)[,1:5])

標準化された行列で特異値分解をした結果とprcomp()関数の出力結果を比較する。

> head(round(Z.svd.score,3)[,1:5])

       [,1]   [,2]   [,3]   [,4]   [,5]

[1,] -0.595  0.260  0.268  0.518 -0.520

[2,]  3.213  0.435  1.089 -1.717  0.079

[3,]  1.798 -1.489 -2.021 -0.433  1.478

[4,] -1.161 -1.124  1.607 -1.371  2.359

[5,]  2.841 -2.078 -0.574 -0.500  0.656

[6,]  2.837 -0.149 -1.400  1.052  0.369

> # 主成分得点の一部を表示する。

> head(round(X.pca$x ,3)[,1:5])

                PC1    PC2    PC3    PC4    PC5

Hokkaido     -0.595  0.260  0.268  0.518 -0.520

Aomori-ken    3.213  0.435  1.089 -1.717  0.079

Iwate-ken     1.798 -1.489 -2.021 -0.433  1.478

Miyagi-ken   -1.161 -1.124  1.607 -1.371  2.359

Akita-ken     2.841 -2.078 -0.574 -0.500  0.656

Yamagata-ken  2.837 -0.149 -1.400  1.052  0.369

このように、データを標準化した結果と同じなり、
変数間の単位の差は、調整されている。

ところで、主成分分析の結果として出力される分散は、
元の行列の分散共分散の特異値に対応するのだった。

また、データを標準化したデータの分散共分散行列は、
相関係数行列に一致するのだった。

# 分散共分散に対応する主成分の分散
round(svd(cov(X))$d, 3)

# 相関係数行列に対応する主成分の分散
round(svd(cor(X))$d, 3)

# prcomp()関数の場合。確か、標準偏差の二乗が分散だった。
round((X.pca$sdev)^2, 3)

以下が、上記の実行結果。


> # 分散共分散に対応する主成分の分散
> round(svd(cov(X))$d, 3)
 [1] 312.518 100.764  65.603  53.541  28.354  22.807  15.126   7.267   6.065   3.257
> 
> # 相関係数行列に対応する主成分の分散
> round(svd(cor(X))$d, 3)
 [1] 3.012 1.507 1.349 1.100 0.872 0.629 0.569 0.468 0.295 0.198
> 
> # prcomp()関数の場合。確か、標準偏差の二乗が分散だった。
> round((X.pca$sdev)^2, 3)
 [1] 3.012 1.507 1.349 1.100 0.872 0.629 0.569 0.468 0.295 0.198

さて、相関係数行列に対応する主成分分析を行った場合、
出てきた分散は、非常に面白い性質を持っている。

# とりあえず、平均を取ってみる。
mean(svd(cor(X))$d)

そして、その実行結果。


> mean(svd(cor(X))$d)
[1] 1

平均は「1」。標準化は、平均が「0」で分散が「1」となるような調整であったので、
当然、分散の平均は「1」となるのであるが、個々の値を見ると「1」以上のものもある。

このことから、主成分分析を行うと、元の分散を個々の再配分していることが解る。

さて、今回は、サンプルデータを用いて、
分散共分散行列と相関係数行列の二つの状況に対応する主成分分析を行ったが、
用いたサンプルデータの場合は、どちらが適切であろうか？

データの性質を考えてみると、全ての変数の値は、
ボランティア活動に費やした「1年当たりの活動従事日数」であった。

つまり、変数間の単位は一定であり、変数間に差はあるものの、
その差の大きさも重要な情報である。

したがって、分散共分散行列の方が適しているように思える。

さて、ここまで二回に分けて、主成分分析の仕組みを解説してきたが、
仕組みを理解し、実行したただけでは、まだ、利用できたと言えない。

重要なことは、結果を「読み取る」こと。

次回からは、可視化の方法も含めて、結果を読み取る方法について考えることにする。