Rでの主成分分析の実行
前々回のエントリで学習した永田・棟近教科書の第9章「主成分分析」にのっている計算例を、自分でRにより実行してみることとする。
前半では、教科書の計算例の実行、後半では、Rのprcomp()関数を使うときに注意しなきゃなと思った点をメモしておく。
永田・棟近教科書の第9章「主成分分析」をRで実行してみる
まず、データの入力。
> # データの入力
>
> 生徒NO <- seq(1, 10, 1)
> 国語 <- c(86,71,42,62,96,39,50,78,51,89)
> 英語 <- c(79,75,43,58,97,33,53,66,44,92)
> 数学 <- c(67,78,39,98,61,45,64,52,76,93)
> 理科 <- c(68,84,44,95,63,50,72,47,72,91)
データを確認のため表示させてみる。
> 成績d <- data.frame(生徒NO=生徒NO, 国語=国語, 英語=英語, 数学=数学, 理科=理科)
> print(成績d)
生徒NO 国語 英語 数学 理科
1 1 86 79 67 68
2 2 71 75 78 84
3 3 42 43 39 44
4 4 62 58 98 95
5 5 96 97 61 63
6 6 39 33 45 50
7 7 50 53 64 72
8 8 78 66 52 47
9 9 51 44 76 72
10 10 89 92 93 91
なお、データをテーブルでみる分には上記のほうがいいけど、データの部分を対象として計算する上でいちいちデータフレームの２列目以降を指定するのはダルいので、行列にしておく。
> 成績 <- cbind(国語, 英語, 数学, 理科)
> print(成績) # これは行列
国語 英語 数学 理科
[1,] 86 79 67 68
[2,] 71 75 78 84
[3,] 42 43 39 44
[4,] 62 58 98 95
[5,] 96 97 61 63
[6,] 39 33 45 50
[7,] 50 53 64 72
[8,] 78 66 52 47
[9,] 51 44 76 72
[10,] 89 92 93 91
>
散布図行列と相関行列をみてみる。
> pairs(成績, pch=16, col="blue") # 散布図行列。点は青にしてみた。
> R <- cor(成績) # 相関行列
> print(R) # 相関行列を出力
国語 英語 数学 理科
国語 1.0000 0.9670 0.3761 0.3113
英語 0.9670 1.0000 0.4146 0.3983
数学 0.3761 0.4146 1.0000 0.9721
理科 0.3113 0.3983 0.9721 1.0000
>
この相関行列から、固有値と固有ベクトルを求めてみる。eigen()関数を使う。
>
> eigen(R) # これで相関行列のスペクトル分解が一瞬で行われる。
$values
[1] 2.720733 1.221800 0.052411 0.005056$vectors
[,1] [,2] [,3] [,4]
[1,] -0.4873 0.5273 -0.4990 0.4853
[2,] -0.5105 0.4740 0.5387 -0.4738
[3,] -0.5083 -0.4807 -0.5041 -0.5063
[4,] -0.4935 -0.5159 0.4547 0.5326
次に、主成分分析を実行してみる。べつに上記の固有値分解を先にやっておかなくても分析は行える。
> result <- prcomp(成績, scale=TRUE)
>
> print(result) # これで、主成分にかかる固有ベクトルが分かる。Standard deviations:
[1] 1.64946 1.10535 0.22894 0.07111Rotation:
PC1 PC2 PC3 PC4
国語 -0.4873 -0.5273 -0.4990 0.4853
英語 -0.5105 -0.4740 0.5387 -0.4738
数学 -0.5083 0.4807 -0.5041 -0.5063
理科 -0.4935 0.5159 0.4547 0.5326
この、Rotationというところは固有ベクトルをタテに束ねた行列であり、主成分負荷量ではない点に注意が必要である。
PC1という列にタテに並んでいる数字は、各科目の得点を合成して第1主成分を表すときの、各変数の重み（係数)を表している。つまりこの表は、「タテに見る」ことが多いものだと気をつけること。
というのも、因子分析だと因子負荷行列を出力するが、あれはだいたいヨコに見るものなので混同しないようにしないといけないわけである。
ちなみに、主成分分析をprincomp()という関数で実行することもできるのだが、こちらは固有ベクトルをタテに束ねた行列が「$loadings」という要素に入っているので、なおさら、因子負荷量（主成分負荷量）のことかと勘違いしやすい。
> result2 <- princomp(成績, cor=TRUE, scores=TRUE)
> result2$loadingsLoadings:
Comp.1 Comp.2 Comp.3 Comp.4
国語 -0.487 0.527 -0.499 0.485
英語 -0.511 0.474 0.539 -0.474
数学 -0.508 -0.481 -0.504 -0.506
理科 -0.493 -0.516 0.455 0.533Comp.1 Comp.2 Comp.3 Comp.4
SS loadings 1.00 1.00 1.00 1.00
Proportion Var 0.25 0.25 0.25 0.25
Cumulative Var 0.25 0.50 0.75 1.00
さてprcomp()の結果に戻る。summaryの中にある寄与率、累積寄与率をみてみる。
> summary(result) # これで、寄与率と累積寄与率もみることができる。
Importance of components:
PC1 PC2 PC3 PC4
Standard deviation 1.65 1.105 0.2289 0.07111
Proportion of Variance 0.68 0.305 0.0131 0.00126
Cumulative Proportion 0.68 0.986 0.9987 1.00000
>
これで、第2主成分までみれば、全体の情報量の98.6%が説明できることがわかる。
つぎに主成分得点を表示してみる。
> result$x # これで主成分得点がみれる。
PC1 PC2 PC3 PC4
[1,] -0.7958 -0.85761 -0.1085806 0.123347
[2,] -1.0732 0.34757 0.2741601 0.043378
[3,] 2.4942 -0.32031 0.1822602 -0.103738
[4,] -1.2840 1.76450 -0.1744056 0.007490
[5,] -1.1646 -1.80252 0.1280482 -0.027363
[6,] 2.4800 0.29770 0.0008193 0.066105
[7,] 0.6428 0.67851 0.2953350 0.041382
[8,] 0.6720 -1.34136 -0.3798562 -0.009874
[9,] 0.5173 1.14860 -0.2661938 -0.050872
[10,] -2.4886 0.08492 0.0484133 -0.089856
主成分得点というのは、各サンプルが、主成分PC1〜PC4に対して持っている値であり、言ってみればPC1〜PC4が座標軸だった場合の、それぞれの軸に関する座標である。（この場合、4次元空間になるので図示できないが。）
で、たとえば「サンプル番号1」の人のPC1〜PC4に対する主成分得点と、各主成分に含まれる「国語」の係数（重み）、つまりさっきの$rotationの値をかけ算して合計すると、サンプル番号1の人の国語の得点を標準化した値になる。
検算のためにやってみる。
> scale(成績)[1,1] # scale()関数はデータを標準化する関数。標準化した成績表の1行1列の値を取ってる。
国語
0.9541
>
> sum(result$x[1,] * result$rotation[1,]) # 上の国語の標準化得点に一致
[1] 0.9541
prcomp()関数を使うときの注意点
さて、Rのprcomp()関数を使うときの注意点として、さしあたり私が学んだことが３つある。１つは既に述べたとおり、主成分の固有ベクトルが出力されている表を主成分負荷量と間違えないことだ。
んで、２つ目は「標準化」をするようにきちんと設定するということで、３つ目は、主成分得点の計算アルゴリズムが他の関数と異なるということである。以下、２つ目と３つ目についてメモっておく。
prcomp()はデフォルトでは、データを標準化しないで分析を行う。つまり、主成分分析には「相関係数行列」から出発する場合と「分散共分散行列」から出発する場合の２種類があると言われる場合の、後者のパターンである。
標準化するのであれば、scaleという引数にscale=TRUEを設定しないといけない。さきほどの計算時には、TRUEを設定しておいた。これのデフォルト値がFALSEになっているので、そのままだと標準化されない。
で、問題は、青木先生の掲示板（このスレッド）でも述べられているように、標準化しなくても良いケースなんて、実際にはほとんど無いということである。尺度（単位）が同じで等分散でなければならないからだ。そんな、実務的に有り得ない仮定を置くぐらいだったら、最初から標準化するのを基本にしておけばいいはずである。
考え方としては、「中心化は必ずする」「標準化は、単位系が同じならしなくてもいいが、その場合は分散の大きな変数の影響が大きくなってしまうので、普通はしたほうがいい」ということになる。（参考追加）
ところが、Rのprcomp()関数の仕様をhelp(prcomp)でみてみると、引数scaleのところには、
scale. a logical value indicating whether the variables should be scaled to have unit variance before the analysis takes place. The default is FALSE for consistency with S, but in general scaling is advisable. Alternatively, a vector of length equal the number of columns of x can be supplied. The value is passed to scale.
つまり、S言語（R言語のもとになった、古い言語）のときにデフォルト値はFALSEになっていたので、それを引きずってますが、基本的にはTRUEにしとくのがオススメですよとわざわざ書いてある。
主成分分析の結果に関する注意点としては、下記のページに関数間の比較として詳しい解説がのってるのだが、ここで３つ目の問題が解説されてる。
湿地福利院|湿地福利院有限公司
分析プログラムにより結果が違う，そんな馬鹿な - 裏 RjpWiki
標準化するしないは単に設定すればすむ問題だが、主成分得点を求める際に、主成分得点の不偏分散が固有値に一致するように計算される（prcomp関数）か、主成分得点の標本分散が固有値に一致するように計算される（princomp関数）かという違いがあるようである。
これはややこしい！
