Rでの主成分分析の結果の違いについて

1+

R には主成分分析のための関数が、
主なものとして、prcomp, princomp,
そして FactoMineR の PCAと、3つあります。

計算結果の数値の違いについて、まとめてみました。

なお、このサイトの主成分分析の説明では、
prcomp の結果に合わせてあります。

Statistics: 主成分分析

ここでは説明の便宜のために、
主成分分析の結果の格納先を、
関数名と同じ名前の、変数とします。

以下を実行済として、説明を進めていきます。

library(FactoMineR);
prcomp<-prcomp(attitude, scale=T);
princomp<-princomp(attitude, cor=T);
pca<-PCA(attitude, scale=T, ncp=7);

なお、データセットは attitude を用い、
また相関行列を用いた計算を行うこととします。

主成分得点の正負については、固有ベクトルに
起因するものであるため、無視して頂いて結構です。

(1)主成分得点の違いについて

princomp$scores と、 pca$ind$coord で
得られる主成分得点の値は同じです。
これらの値と、 prcomp$x で得られる
主成分得点の値は、微妙に異なります。

cbind(cbind(prcomp$x[,1], princomp$scores[,1]), pca$ind$coord[,1]);

これは、prcomp で得られる主成分得点が、
「主成分得点の不偏分散 = 固有値」
となっているのに対して、
princomp と pca で得られる主成分得点が、
「主成分得点の標本分散 = 固有値」
となっていることに起因する違いです。

matrix(c(eigen(cor(attitude))$values, diag(var(prcomp$x)), diag(var(princomp$scores)), diag(var(pca$ind$coord))), 7, 4);

具体的に、どこで不偏/標本分散を使っているのかというと。。。
主成分得点算出の以下の式の、

s<-scale(attitude) %*% t(solve(eigen(cor(attitude))$vectors));

ここの、 scale関数、
もとデータの標準化の部分で、違いが出ています。

(ここでは、相関行列をもとにしているため、
上式後半の、固有ベクトルには、標本/不偏分散は影響ありません。)

scale関数の標準化の際に、標準偏差=分散の平方根で
割っているため、この分散が不偏分散か、標本分散か、で
違いが出ています。

そのため、princomp, pca と同じ値を求めるためには、
以下のように、標本分散の平方根を用いた標準化を
行う必要があります。

(scale(attitude)/sqrt((nrow(attitude)-1)/nrow(attitude))) %*% t(solve(eigen(cor(attitude))$vectors));

scale 関数は、標準化された値、すなわち
各列の偏差を、各列の標準偏差で割った値を出力します。

scale(attitude)[,1];
# 標準化の scale 関数の中身は以下
(attitude[,1]-mean(attitude[,1]))/sqrt(var(attitude[,1]));
# 不偏分散の var 関数を標本分散に変換
(v<-(nrow(attitude)-1)/nrow(attitude));
(attitude[,1]-mean(attitude[,1]))/sqrt(var(attitude[,1])*v);
# 以下式変形。 scale を v の平方根で割ることになる
(attitude[,1]-mean(attitude[,1]))/(sqrt(var(attitude[,1]))*sqrt(v));
(attitude[,1]-mean(attitude[,1]))/(sqrt(var(attitude[,1])))/sqrt(v);
(scale(attitude)/sqrt(v))[,1];

そのため、prcomp の主成分得点を、
(N-1 / N)の平方根で割れば、
princomp と pcaの結果と一致します。

逆もしかり、
princomp と pca の主成分得点を、
(N / N-1)の平方根で割れば、
prcomp の結果と一致します。
(Nはデータ件数を指します)

(v<-(nrow(attitude)-1)/nrow(attitude));
cbind(cbind(prcomp$x[,1]/sqrt(v), princomp$scores[,1]), pca$ind$coord[,1]);
cbind(cbind(prcomp$x[,1], princomp$scores[,1]/sqrt(1/v)), pca$ind$coord[,1]/sqrt(1/v));

(2)主成分負荷量について

FactoMineR のPCAで矢印プロットされる、主成分負荷量。
これは、prcomp や princomp の結果の
オブジェクトの中には含まれません。

自分で計算して、値を求める必要があります。

主成分負荷量は、因子分析の因子負荷量と同じく、
固有値の平方根に、固有ベクトルをかけたもの。

これで Y先生いわく、軽みづけ、
重み付けでなく、固有値の低いものを小さくするための
ウェートをかけているのだそうです。

cbind((prcomp$rotation %*% diag(prcomp$sdev))[,1], pca$var$coord[,1]);
cbind((princomp$loadings %*% diag(princomp$sdev))[,1], pca$var$coord[,1]);
cbind((eigen(cor(attitude))$vectors %*% diag(sqrt(eigen(cor(attitude))$values)))[,1], pca$var$coord[,1]);

prcomp$rotation また princomp$loadings は
固有ベクトルそのものを格納しています。
sdev は固有値の平方根を格納しています。
また、3行目は固有値から計算しています。
全て、同じ結果を返します。

以上です。


参考資料

http://aoki2.si.gunma-u.ac.jp/R/princomp2.html
分析プログラムにより結果が違う,そんな馬鹿な

カテゴリー:


1+