StatsBeginner: 初学者の統計学習ノート

初学者が統計学、機械学習、R、Pythonの勉強の過程をメモっていくノート。

Rのprcomp()関数で主成分分析をするときの注意点

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

 なお、データをテーブルでみる分には上記のほうがいいけど、データの部分を対象として計算する上でいちいちデータフレームの2列目以降を指定するのはダルいので、行列にしておく。

> 成績 <- 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") # 散布図行列。点は青にしてみた。

f:id:midnightseminar:20140727113807p:plain:w400

> 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.07111

Rotation:
    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$loadings

Loadings:
   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.533

         Comp.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()関数を使うときの注意点として、さしあたり私が学んだことが3つある。1つは既に述べたとおり、主成分の固有ベクトルが出力されている表を主成分負荷量と間違えないことだ。
 んで、2つ目は「標準化」をするようにきちんと設定するということで、3つ目は、主成分得点の計算アルゴリズムが他の関数と異なるということである。以下、2つ目と3つ目についてメモっておく。


 prcomp()はデフォルトでは、データを標準化しないで分析を行う。つまり、主成分分析には「相関係数行列」から出発する場合と「分散共分散行列」から出発する場合の2種類があると言われる場合の、後者のパターンである。
 標準化するのであれば、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にしとくのがオススメですよとわざわざ書いてある。


 主成分分析の結果に関する注意点としては、下記のページに関数間の比較として詳しい解説がのってるのだが、ここで3つ目の問題が解説されてる。


 Statistics: Rでの主成分分析の結果の違いについて
 分析プログラムにより結果が違う,そんな馬鹿な - 裏 RjpWiki


 標準化するしないは単に設定すればすむ問題だが、主成分得点を求める際に、主成分得点の不偏分散が固有値に一致するように計算される(prcomp関数)か、主成分得点の標本分散が固有値に一致するように計算される(princomp関数)かという違いがあるようである。
 これはややこしい!