前回は、特異値分解を通して主成分分析の仕組みについて整理した。
主成分分析は、「対象の本質」と「変数の本質」に分離し、
「対象のみに焦点」を当てて分析する方法だった。
そして、対象の本質のみを記述した部分を主成分得点と呼び、
変数の本質のみを記述した部分を主成分負荷量と呼んだ。
また、対象のみの状態にすることで、変数間の関係から開放され、
分散のみに注目するだけで良い状況になった。
前回の話では、特異値分解を用いて主成分分析を行ったが、
もちろん、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
> 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
主成分分析は、「対象の本質」と「変数の本質」に分離し、
「対象のみに焦点」を当てて分析する方法だった。
そして、対象の本質のみを記述した部分を主成分得点と呼び、
変数の本質のみを記述した部分を主成分負荷量と呼んだ。
また、対象のみの状態にすることで、変数間の関係から開放され、
分散のみに注目するだけで良い状況になった。
前回の話では、特異値分解を用いて主成分分析を行ったが、
もちろん、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=","))
まずは、分散共分散行列、相関係数行列、散布図行列を確認する。
毎回、これらを確認することを忘れてはいけない。
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 というのは、元のデータの真ん中。列平均値。特異値分解とは無関係。
> 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(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)
このように、データを標準化した結果と同じなり、
変数間の単位の差は、調整されている。
ところで、主成分分析の結果として出力される分散は、
元の行列の分散共分散の特異値に対応するのだった。
また、データを標準化したデータの分散共分散行列は、
相関係数行列に一致するのだった。
# 分散共分散に対応する主成分の分散
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年当たりの活動従事日数」であった。
つまり、変数間の単位は一定であり、変数間に差はあるものの、
その差の大きさも重要な情報である。
したがって、分散共分散行列の方が適しているように思える。
さて、ここまで二回に分けて、主成分分析の仕組みを解説してきたが、
仕組みを理解し、実行したただけでは、まだ、利用できたと言えない。
重要なことは、結果を「読み取る」こと。
次回からは、可視化の方法も含めて、結果を読み取る方法について考えることにする。
[,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年当たりの活動従事日数」であった。
つまり、変数間の単位は一定であり、変数間に差はあるものの、
その差の大きさも重要な情報である。
したがって、分散共分散行列の方が適しているように思える。
さて、ここまで二回に分けて、主成分分析の仕組みを解説してきたが、
仕組みを理解し、実行したただけでは、まだ、利用できたと言えない。
重要なことは、結果を「読み取る」こと。
次回からは、可視化の方法も含めて、結果を読み取る方法について考えることにする。
0 件のコメント:
コメントを投稿