「因子軸を回転って結局どういうイメージですか」
という質問をされたのですね。
探索的因子分析の後は,解釈をしやすくするために軸を回転しますよ,と説明しているのだけど,そのイメージがわかりにくいということなので,Rでプロットしてみることにしました。
psychパッケージのデータセットbfiを例にあげます。これは五次元のデータなんですけど,仮に2因子だとしましょう。
(後で出てきますが3次元空間の回転行列から角度を求めるのは行列計算的にちょっと面倒なのです。)
因子数を2にし,回転なし,バリマックス回転,プロマックス回転の結果をオブジェクトに代入しておきます。
library(psych)
library(ggplot2)
## ## 次のパッケージを付け加えます: 'ggplot2' ## ## 以下のオブジェクトは 'package:psych' からマスクされています: ## ## %+%
dat <- bfi[1:25]
fa.none <- fa(dat,2,rotate="none")
fa.vari <- fa(dat,2,rotate="varimax")
fa.prom <- fa(dat,2,rotate="promax")
直交回転の場合
さてここから加工。まず因子負荷行列をmatrix型で抜き出し,データフレーム型に変換。psychパッケージのfaオブジェクトはいろいろ付いているからね。無駄なものをそぎ落とします。
none.ld <- matrix(fa.none$loadings,ncol=fa.none$factors)
rotated <- matrix(fa.vari$loadings,ncol=fa.none$factors)
none.ld.df <- as.data.frame(none.ld)
rotated.df <- as.data.frame(rotated)
因子名をわかりやすくFx,Fyにしましょう(2次元だから)
colnames(none.ld.df) <- c("Fx","Fy")
colnames(rotated.df) <- c("Fx","Fy")
回転なしの因子負荷量をつかって因子負荷行列のプロットを描いてみます。
g <- ggplot(none.ld.df,aes(x=Fx,y=Fy))+geom_point() +xlim(-1,1)+ylim(-1,1)
g <- g + geom_hline(yintercept=0,colour="blue") + geom_vline(xintercept=0,colour="blue")
g
回転後の因子負荷行列の散布図も描いてみます。
g <- ggplot(rotated.df,aes(x=Fx,y=Fy))+geom_point() +xlim(-1,1)+ylim(-1,1)
g <- g + geom_hline(yintercept=0,colour="blue") + geom_vline(xintercept=0,colour="blue")
g
うむ,確かにこれでは分かりにくい。
さて,回転は行列をかけることで行っています。座標を変換するんですね。Wikiのページにもあるように,回転行列とはcos,sinをつかった行列です。つ回転行列
この回転行列はrot.matの中に入っています。
rm <- fa.vari$rot.mat
1行1列目は$\cos \theta$の値でしたね。この値から逆に角度を求めます。逆cosはacosです。
rad <- acos(rm[1,1])
rad
## [1] 0.7139151
0.71ラジアンぐらいですね
図にしましょう。まず回転前のプロット。
g <- ggplot(none.ld.df,aes(x=Fx,y=Fy))+geom_point() +xlim(-1,1)+ylim(-1,1)
g <- g + geom_hline(yintercept=0,colour="blue") + geom_vline(xintercept=0,colour="blue")
回転軸を引きます。角度を傾きにするため,tanをとっています
g <- g + geom_abline(intercept=0,slope=tan(rad),colour="red",linetype=2)
第2次元は第1次元と直角に交わるので,0.5piを足してtanをとっています。
g <- g + geom_abline(intercept=0,slope=tan(rad+(pi/2)),colour="red",linetype=2)
# 描画
g
たしかに,元の青い軸のプロットよりも赤い軸のプロットの方が,点と軸の距離が近いようです。強調した,というのがちょっとは実感できるでしょうか。
斜交回転
promax回転の場合。最初の処理は同じです。
none.ld <- matrix(fa.none$loadings,ncol=fa.none$factors)
rotated <- matrix(fa.prom$loadings,ncol=fa.none$factors)
none.ld.df <- as.data.frame(none.ld)
rotated.df <- as.data.frame(rotated)
colnames(none.ld.df) <- c("Fx","Fy")
colnames(rotated.df) <- c("Fx","Fy")
回転行列を抜き出します。
rm <- fa.prom$rot.mat
角度を抜き出します。
rad <- acos(rm[1,1])
さて,斜交回転の場合は軸が相関します。相関行列はPhiに入っています。相関係数は(\cos \theta) の結果なので,逆算(acos)して第一軸と第二軸の角度を求めます。
phi_rad <- acos(fa.prom$Phi[1,2])
これでプロット。
# 回転前の散布図
g <- ggplot(none.ld.df,aes(x=Fx,y=Fy))+geom_point() +xlim(-1,1)+ylim(-1,1)
g <- g + geom_hline(yintercept=0,colour="blue") + geom_vline(xintercept=0,colour="blue")
# 第1次元を引きます。角度をtanしたのが傾斜です
g <- g + geom_abline(intercept=0,slope=tan(rad),colour="red",linetype=2)
# 第2次元は第1次元に角度分たしたのが傾斜になります。
g <- g + geom_abline(intercept=0,slope=tan(rad+phi_rad),colour="red",linetype=2)
# 描画
g
こんなもんでどうでしょう。
本当のこというと,適当な2軸を選んで描画し,軸をぐいっと動かして相関係数を見せるようなインタラクティブなR画面を作りたかったけど,それは技量と時間の限界ということで,誰か助けてエロい人。