R には主成分分析のための関数が、
主なものとして、prcomp, princomp,
そして FactoMineR の PCAと、３つあります。
計算結果の数値の違いについて、まとめてみました。
なお、このサイトの主成分分析の説明では、
prcomp の結果に合わせてあります。
ここでは説明の便宜のために、
主成分分析の結果の格納先を、
関数名と同じ名前の、変数とします。
以下を実行済として、説明を進めていきます。
library(FactoMineR);
prcomp<-prcomp(attitude, scale=T);
princomp<-princomp(attitude, cor=T);
pca<-PCA(attitude, scale=T, ncp=7);
なお、データセットは attitude を用い、
また相関行列を用いた計算を行うこととします。
主成分得点の正負については、固有ベクトルに
起因するものであるため、無視して頂いて結構です。
（１）主成分得点の違いについて
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));
（２）主成分負荷量について
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 は固有値の平方根を格納しています。
また、３行目は固有値から計算しています。
全て、同じ結果を返します。
以上です。
参考資料
http://aoki2.si.gunma-u.ac.jp/R/princomp2.html
分析プログラムにより結果が違う，そんな馬鹿な