因子分析におけるMAP基準について

はじめに

探索的因子分析において,因子数を決める際にはさまざまな基準があるのですが,その中にMAP基準というものがあります。必要があって調べた際のメモを残しておきます。

MAP基準とは

MAPはMinimum average partical correlationの略です。オリジナルはVelicer(1976)による次の論文です。https://doi-org/10.1007/BF02293557

$p$個の観測変数の相関行列を$R$とします。サイズは$p \times p$です。主成分分析の結果である$p \times m$のパターン行列を$A$とします。そのとき,偏分散共分散行列$C^*$は

$$C^* = R- AA^\prime$$

と書けます。ここからさらに偏相関係数の行列は

$$ R^* = D^{-1/2}( R- AA^\prime)D^{-1/2} \tag{1} $$

と書けます。ここで $D$は,

$$D = Diag(C^*) = Diag(R- AA^\prime)$$

です。ここで,$R^*$のi行目,j列目の要素を $r ^* _{ij}$としたときに,偏相関行列の2乗の値の平均が

$$ f_m = \frac{\sum \sum_{i \neq j} (r ^* _{ij})^2 }{p(p-1)} \tag{2}$$

として表されます。MAP基準とは,パターン行列のmの値を変えていきながらこの値を計算して,最小になる$f_m$の$m$を因子数として採択しようという方法です。

実際の数値例を確認

実際にデータを触ってみながら,手順を確認してみます。psychvss関数の中身を覗きながら書きました。使うデータはpsychパッケージの中に入っているThurstone.9という9つのテストの相関行列です。

library(psych)

data("Thurstone.9")
R <- Thurstone.9
Prefixes Suffixes Vocabulary Sentences First.Last FirstLetters FourLetters Completion SameorOpposite
Prefixes 1.000 0.554 0.227 0.189 0.461 0.506 0.408 0.280 0.241
Suffixes 0.554 1.000 0.296 0.219 0.479 0.530 0.425 0.311 0.311
Vocabulary 0.227 0.296 1.000 0.769 0.237 0.243 0.304 0.718 0.730
Sentences 0.189 0.219 0.769 1.000 0.212 0.226 0.291 0.681 0.661
First.Last 0.461 0.479 0.237 0.212 1.000 0.520 0.514 0.313 0.245
FirstLetters 0.506 0.530 0.243 0.226 0.520 1.000 0.473 0.348 0.290
FourLetters 0.408 0.425 0.304 0.291 0.514 0.473 1.000 0.374 0.306
Completion 0.280 0.311 0.718 0.681 0.313 0.348 0.374 1.000 0.672
SameorOpposite 0.241 0.311 0.730 0.661 0.245 0.290 0.306 0.672 1.000

とりあえず,因子数が1のときの$f_m$から求めてみましょう。 まず,この相関行列を固有値分解して固有値ベクトルを求め,そこからパターン行列をえます。

e <- eigen(R)
evect <- e$vectors
comp <- evect %*% diag(sqrt(e$values))
round(comp, 3)
       [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]
 [1,] -0.595  0.489 -0.385  0.308  0.067  0.396  0.042  0.015  0.002
 [2,] -0.645  0.441 -0.356  0.057  0.014 -0.495 -0.102 -0.053 -0.043
 [3,] -0.763 -0.505 -0.054  0.036  0.057 -0.037 -0.076  0.088  0.375
 [4,] -0.716 -0.526  0.022  0.025  0.037  0.072 -0.317  0.222 -0.230
 [5,] -0.618  0.471  0.312 -0.193  0.510  0.016  0.030  0.022 -0.001
 [6,] -0.646  0.463 -0.022 -0.461 -0.358  0.132 -0.039  0.083  0.037
 [7,] -0.645  0.326  0.522  0.362 -0.258 -0.077  0.025  0.021  0.008
 [8,] -0.787 -0.369  0.041 -0.063 -0.027  0.094 -0.035 -0.475 -0.054
 [9,] -0.748 -0.433 -0.086 -0.043 -0.020 -0.081  0.458  0.123 -0.108

続いて,(1)式から偏相関係数行列を求めます。$D$は対角行列ですので,$D^{-1/2}$は$D$の対角成分の平方根の逆数を対角成分に持つ対角行列です。

Cstar <- R - comp[,1] %*% t(comp[,1])
d <- diag(1/sqrt(diag(Cstar)))
Rstar <- d %*% Cstar %*% d
round(Rstar,3)

       [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]
 [1,]  1.000  0.277 -0.438 -0.424  0.148  0.198  0.039 -0.380 -0.384
 [2,]  0.277  1.000 -0.397 -0.456  0.134  0.194  0.015 -0.416 -0.338
 [3,] -0.438 -0.397  1.000  0.493 -0.461 -0.507 -0.381  0.295  0.371
 [4,] -0.424 -0.456  0.493  1.000 -0.420 -0.445 -0.321  0.272  0.270
 [5,]  0.148  0.134 -0.461 -0.420  1.000  0.201  0.192 -0.356 -0.416
 [6,]  0.198  0.194 -0.507 -0.445  0.201  1.000  0.096 -0.341 -0.383
 [7,]  0.039  0.015 -0.381 -0.321  0.192  0.096  1.000 -0.284 -0.349
 [8,] -0.380 -0.416  0.295  0.272 -0.356 -0.341 -0.284  1.000  0.203
 [9,] -0.384 -0.338  0.371  0.270 -0.416 -0.383 -0.349  0.203  1.000

(2)式に従って,この行列の非対角要素の2乗の平均をとります。

p <- ncol(R)
diag(Rstar) <- 0
sum(Rstar^2)/(p*(p-1))
[1] 0.1142189

ということで,因子数が1のときの$f_m$が求まりました。

あとは,この手続きを因子数を増やしながら繰り返していくだけなのですが,一応因子数が2の場合もやってみましょう。といってもパターン行列は求まっているのでそこの部分を変えて計算するだけです。

Cstar2 <- R - comp[,1:2] %*% t(comp[,1:2])
d2 <- diag(1/sqrt(diag(Cstar2)))
Rstar2 <- d2 %*% Cstar2 %*% d2
diag(Rstar2) <- 0
sum(Rstar2^2)/(p*(p-1))
[1] 0.03953237

さて,$f_1=0.114$で$f_2=0.040$ですので,比較すると$m=2$の方が小さい値です。これをさらに先まで繰り返して,一番小さい値を出す因子数を採択することになります。

psychの関数

psychvss関数は,この手続きを全て自動でやってくれます。デフォルトでは8因子までですが,引数にn=10などと指定すれば,もっと多い因子数までの値を計算してくれます(あまり必要になることはないでしょうが)。

res.vss <- vss(R)
res.vss$map
[1] 0.11421886 0.03953237 0.06829390 0.12575972 0.20558816 0.35901589 0.54588031 1.00000000

この中だと3つめの値が一番小さいので,3因子を選択することとなります。