Rの主成分分析パッケージ「principal」と「prcomp」の違いをまとめてみた。

Rの主成分分析パッケージ「principal」と「prcomp」の違いをまとめてみた。
Rで主成分分析をする際には、psychパッケージのprincipalとRデフォルトのprcompを使うことが多いと思います。しかし、この2つは微妙に違う結果を出してきます。そして、その内容をいつも忘れてしまう。。。というわけで、今回は備忘録を兼ねて、2つの違い説明していきます。ちなみに、同様の記事が旧サイトにありますが、決してパクリではないので。

準備

サンプルデータの作成

先ず、サンプルデータを作成します。主成分分析は相関が高い変数をまとめて、新しい合成変数(主成分)を作成する手法なので、今回はmvrnormを利用して、相関がある疑似データを作成してみます。

In [1]:
#桁数の設定
options(digits=3) 

#MASSパッケージの呼び出し
library(MASS)

#擬似的に作成するデータセットの相関を設定
set.seed(100)
dammy_cor = matrix(c(1.0, 0.8, 0.8, 0.3, 0.3, 0.3, 0.3,
                     0.8, 1.0, 0.8, 0.3, 0.3, 0.3, 0.3,
                     0.8, 0.8, 1.0, 0.3, 0.3, 0.3, 0.3,
                     0.3, 0.3, 0.3, 1.0, 0.8, 0.8, 0.8,
                     0.3, 0.3, 0.3, 0.8, 1.0, 0.8, 0.8,
                     0.3, 0.3, 0.3, 0.8, 0.8, 1.0, 0.8,
                     0.3, 0.3, 0.3, 0.8, 0.8, 0.8, 1.0), ncol=7)

#疑似データの作成
samp_dat <- mvrnorm(n=500, mu=c(0,0,0,0,0,0,0), Sigma=dammy_cor, empirical=F)
colnames(samp_dat) <- paste('X', 1:7, sep="") #ヘッダー付与
cor(samp_dat)
X1 X2 X3 X4 X5 X6 X7
X1 1.000 0.810 0.802 0.264 0.255 0.281 0.245
X2 0.810 1.000 0.835 0.274 0.260 0.285 0.255
X3 0.802 0.835 1.000 0.280 0.243 0.281 0.249
X4 0.264 0.274 0.280 1.000 0.774 0.815 0.771
X5 0.255 0.260 0.243 0.774 1.000 0.802 0.800
X6 0.281 0.285 0.281 0.815 0.802 1.000 0.823
X7 0.245 0.255 0.249 0.771 0.800 0.823 1.000

主成分分析の実行

In [2]:
principal_result <- psych::principal(samp_dat,nfactors=7,rotate='none')
prcomp_result <- prcomp(samp_dat,scale=TRUE)

それぞれのパッケージを使って主成分分析を実行して、結果のオブジェクトを格納しておきます。

主成分負荷量

主成分負荷量は、主成分と各変数との相関であり、主成分の解釈する際に用います。

principalの主成分負荷量

principalの場合は、結果のloadingsで確認することができます。

In [3]:
round(principal_result$loadings[],5)
PC1 PC2 PC3 PC4 PC5 PC6 PC7
X1 0.635 0.678 0.0916 0.2214 -0.27678 0.0550 0.00612
X2 0.647 0.684 0.0142 -0.0746 0.14814 -0.1845 0.22697
X3 0.641 0.686 -0.0910 -0.1464 0.12911 0.1375 -0.23031
X4 0.829 -0.378 -0.3694 0.1032 0.00732 0.0993 0.11194
X5 0.823 -0.405 0.2212 0.2284 0.23419 0.0207 -0.05585
X6 0.851 -0.385 -0.0468 -0.0709 -0.13234 -0.2727 -0.16753
X7 0.825 -0.411 0.1881 -0.2564 -0.10774 0.1565 0.11236

procompの主成分負荷量

procompの場合、分析結果として主成分負荷量は出力されません。別途、計算が必要です。

In [4]:
round(t(t(prcomp_result$rotation)*prcomp_result$sdev)[,drop = FALSE],5)
PC1 PC2 PC3 PC4 PC5 PC6 PC7
X1 -0.635 0.678 0.0916 -0.2214 0.27678 0.0550 -0.00612
X2 -0.647 0.684 0.0142 0.0746 -0.14814 -0.1845 -0.22697
X3 -0.641 0.686 -0.0910 0.1464 -0.12911 0.1375 0.23031
X4 -0.829 -0.378 -0.3694 -0.1032 -0.00732 0.0993 -0.11194
X5 -0.823 -0.405 0.2212 -0.2284 -0.23419 0.0207 0.05585
X6 -0.851 -0.385 -0.0468 0.0709 0.13234 -0.2727 0.16753
X7 -0.825 -0.411 0.1881 0.2564 0.10774 0.1565 -0.11236

ちなみに、結果で出力されてくるrotationは固有ベクトルです。結果が似ているので、主成分負荷量と間違えないように注意して下さい。固有ベクトルは主成分得点を計算する際の重みになります。

固有値

固有値は各主成分の分散になります。よく固有値が1以上(ガットマン基準)というのは、標準化した変数は分散が1になるので、それより分散が大きければデータを分類する軸として意味あるよねって解釈です。

principalの固有値

principalの場合は、valuesで結果を返してくれます。

In [5]:
round(principal_result$values,5)
  1. 4.00399
  2. 2.02123
  3. 0.23986
  4. 0.20965
  5. 0.19924
  6. 0.1651
  7. 0.16094

procompの固有値

procompの場合は、固有値を返してくれません。主成分の標準偏差の結果sdevを使って、別途計算します。標準偏差の自乗が分散なので。

In [6]:
round(prcomp_result$sdev^2,5)
  1. 4.00399
  2. 2.02123
  3. 0.23986
  4. 0.20965
  5. 0.19924
  6. 0.1651
  7. 0.16094

ちゃんと固有値の値は同じなりますね。ちなみに、当然ですが、eigen関数を使って固有値を確認してみても同じになります。

In [7]:
round(eigen(cor(samp_dat))$values,5)
  1. 4.00399
  2. 2.02123
  3. 0.23986
  4. 0.20965
  5. 0.19924
  6. 0.1651
  7. 0.16094

主成分得点

固有ベクトルと各変数の値を掛け合わせることで、主成分上におけるスコアを計算することができます。

principalの主成分得点

scoresというかたちで、主成分得点の結果を出力してくれます。

In [8]:
head(principal_result$scores) #上位5レコードの主成分得点
PC1 PC2 PC3 PC4 PC5 PC6 PC7
0.4503 1.462 -1.148 0.534 -1.237 -1.047 0.607
-0.1680 -0.244 1.943 -0.151 -1.492 -0.221 -0.272
0.0374 0.386 0.672 1.402 0.948 0.767 0.632
-0.9372 1.881 -0.274 1.989 -0.655 1.289 0.289
-0.1536 -0.157 0.295 0.793 -0.542 -1.784 -0.226
-0.3611 0.399 1.378 1.201 -0.660 0.350 -0.995

procompの主成分得点

xが主成分得点になります。

In [9]:
head(prcomp_result$x)
PC1 PC2 PC3 PC4 PC5 PC6 PC7
-0.9010 2.078 -0.562 -0.244 0.552 -0.4254 -0.2434
0.3362 -0.347 0.951 0.069 0.666 -0.0899 0.1090
-0.0749 0.548 0.329 -0.642 -0.423 0.3116 -0.2534
1.8753 2.674 -0.134 -0.911 0.292 0.5236 -0.1159
0.3073 -0.224 0.145 -0.363 0.242 -0.7248 0.0909
0.7225 0.567 0.675 -0.550 0.294 0.1424 0.3990

スコアが異なってます。最初、見たときになんでって思いましたよw実は、principalの方は出力された主成分得点を標準化したものを結果として返してきます。なので、prcompの結果と一致はしません。prcompの結果を標準化してみます。

In [10]:
head(scale(prcomp_result$x))
PC1 PC2 PC3 PC4 PC5 PC6 PC7
-0.4503 1.462 -1.148 -0.534 1.237 -1.047 -0.607
0.1680 -0.244 1.943 0.151 1.492 -0.221 0.272
-0.0374 0.386 0.672 -1.402 -0.948 0.767 -0.632
0.9372 1.881 -0.274 -1.989 0.655 1.289 -0.289
0.1536 -0.157 0.295 -0.793 0.542 -1.784 0.226
0.3611 0.399 1.378 -1.201 0.660 0.350 0.995

符号の違いはありますが、無事に一致しました。

軸の回転¶

principalは、GPArotationを呼び出すことで、軸の回転を行うことができます。回転を行うことで、軸の解釈がしやすくなります。統計的に良いかという議論はひとまず置いておきます。

In [11]:
library(psych)
library(GPArotation)

#主成分分析の結果に対してバリマックス回転を行う
principal_result_rotate <- psych::principal(samp_dat,nfactors=7,rotate='varimax')

回転前の主成分負荷量

In [12]:
#主成分負荷量
principal_result_rotate$loadings[]
RC2 RC3 RC4 RC1 RC6 RC5 RC7
X1 0.866 0.0809 0.0721 0.0844 0.0582 0.4766 -0.003033
X2 0.944 0.0892 0.0830 0.0873 0.0579 -0.0687 -0.280488
X3 0.940 0.1070 0.0795 0.0617 0.0562 -0.0730 0.292586
X4 0.142 0.8570 0.3295 0.3242 0.1770 0.0223 0.003079
X5 0.120 0.3953 0.3716 0.8094 0.1882 0.0275 -0.002570
X6 0.146 0.4803 0.4351 0.3855 0.6396 0.0301 0.000605
X7 0.118 0.3895 0.8165 0.3602 0.1936 0.0221 0.000372

回転後の主成分負荷量

In [13]:
#主成分得点(回転後)
head(principal_result_rotate$scores)
RC2 RC3 RC4 RC1 RC6 RC5 RC7
1.3342 0.6916 -0.856 -1.213 0.7311 1.018 -0.951
-0.4470 -1.5622 1.127 0.139 0.6323 1.374 -0.123
0.3136 -0.0957 -0.546 1.254 -1.4778 0.410 -0.347
0.8272 -0.1929 -1.644 -0.508 -1.7786 1.917 0.319
-0.3050 -0.4142 -0.805 0.301 1.4316 0.615 -0.928
-0.0434 -1.2448 -0.387 0.789 -0.0473 1.538 0.697

回転後の主成分得点は、回転前の主成分得点にrot.matを掛けた値になります。

In [14]:
head(principal_result$scores %*% principal_result_rotate$rot.mat)
-1.213 1.3342 -0.6916 0.856 -1.018 -0.7311 0.951
0.139 -0.4470 1.5622 -1.127 -1.374 -0.6323 0.123
1.254 0.3136 0.0957 0.546 -0.410 1.4778 0.347
-0.508 0.8272 0.1929 1.644 -1.917 1.7786 -0.319
0.301 -0.3050 0.4142 0.805 -0.615 -1.4316 0.928
0.789 -0.0434 1.2448 0.387 -1.538 0.0473 -0.697

軸の順序、符号がズレてますが、値は一致してます。

まとめ

psychパッケージのprincipalと標準のprcompの結果の比較してみました。個人的にはprincipalの方が、固有値や主成分負荷量を出力してくれるので、便利な印象です。あと、禁じ手?の軸の回転もできますし。prcompは固有ベクトルと主成分負荷量と出力の形が同じで、間違えてしまいそう。利用する場合はご注意を!

 

Rカテゴリの最新記事