あざらしの日常

あざらしが腹を叩く片手間で書いたメモ

因子分析の情報量規準

とあるところで因子分析の情報量規準について見て、少し思うことがあったので簡単にまとめた。

因子分析のモデルは一般的には特異モデルであるため、BICの仮定が成立していない可能性がある。
（因子分析におけるBICの妥当性については後で述べるとして…）こういった場合、特異モデルの情報量規準として

  1. WBIC (Watanabe, 2013)
  2. sBIC (Drton & Plummer, 2017)

を使うことを検討することになるであろう。

因子分析に用いる場合、この2つでどちらの方が良さそうだろうか？
この2つの比較であれば以下の観点でsBICに軍配が上がると思われる。

  • 対数周辺尤度の漸近展開の主要2項からのバイアス

因子分析であれば実対数閾値の理論値が知られている。したがってsBICの方が対数周辺尤度の漸近展開の主要2項はよく近似されると期待される。一方で、WBICはサンプルサイズが小さいときには対数周辺尤度の漸近展開の主要2項からのバイアスが小さくないことが知られている。
sBICに不利な点があるとするならば最適なパラメータのところをMLEで置き換えているため、それによるバイアスが入りうるところだと考えられる。

  • ばらつき

WBICはMCMCにより計算されるためばらつきが出てしまう。一方でsBICであればそのようなばらつきは出ない。

  • 計算量

よく知られている通りMCMCの計算には非常に時間がかかってしまうためWBICの計算も大変になってしまう。sBICであればすぐに計算が可能である。

また、因子分析であればRのsBICパッケージがあるので簡単に計算できる。
実際に簡単なケースでシミュレーションをしてみよう。
シミュレーションの設定は、
Fitting a Bayesian Factor Analysis Model in Stan
を少し変えたものにしておく。具体的には無理やり特異的なケースにするため、
l1 <- c(0.99, 0.00, 0.25, 0.00, 0.80, 0.00, 0.50, 0.00, 0.00, 0.00)
l2 <- c(0.00, 0.90, 0.25, 0.40, 0.00, 0.50, 0.00, 0.00, -0.30, -0.30)
l3<- c(0.00, 0.00, 0.85, 0.80, 0.00, 0.75, 0.75, 0.00, 0.80, 0.80)
のl3をl2に似せて
l3 <- c(0.00, 0.95, 0.25, 0.40, 0.00, 0.50, 0.00, 0.00, -0.30, -0.35)
としておこう。こうすることで、実際は潜在因子が3つであるが2つに近い状態に見えてしまうであろう。
（因子分析的にこの設定に意味があるとは思えないが、元の設定からアドホックに特異なケースにしたものです）
簡単に設定をまとめておくと、真の因子は3つ（そのうちの2つは似てる）で、観測されている変数は10、サンプルサイズを50とし、因子数が1～5のうちどれが最も妥当であるかを周辺尤度を基準に決めようというもの。

WBICの計算は上のリンクのStanコードをちょっといじって計算させました（が、リンク元のコード部分が（ちょこちょこおかしなところがあり）正しい保証がないので、WBICの計算部分は今回は載せないことにしました）
後から気づきましたが、因子分析のstanコードは以下を参考にするのが良さそうです。
StanなMCMCで探索的因子分析やってみた | Kosugitti Labo ver.9

library(MASS)
library(sBIC)

D <- 3
P <- 10 
N <- 50

mu_theta <- rep(0,D) 
mu_epsilon <- rep(0,P) 
Phi <- diag(rep(1,D))
Psi <- diag(c(0.2079, 0.19, 0.1525, 0.20, 0.36, 0.1875, 0.1875, 1.00, 0.27, 0.27))
l1 <- c(0.99, 0.00, 0.25, 0.00, 0.80, 0.00, 0.50, 0.00, 0.00, 0.00)
l2 <- c(0.00, 0.90, 0.25, 0.40, 0.00, 0.50, 0.00, 0.00, -0.30, -0.30)
l3 <- c(0.00, 0.95, 0.25, 0.40, 0.00, 0.50, 0.00, 0.00, -0.30, -0.35)
L <- cbind(l1,l2,l3) 

fa <- FactorAnalyses(numCovariates = 10, maxNumFactors = 5)

iter <- 100

res <- data.frame(model_BIC=numeric(iter), model_sBIC=numeric(iter), model_WBIC=numeric(iter))

for(i in 1:iter){
  Theta <- mvrnorm(N, mu_theta, Phi) 
  Epsilon <- mvrnorm(N, mu_epsilon, Psi) 
  Y <- Theta%*%t(L) + Epsilon
  
  res.fa <- sBIC(Y, fa)
  res$model_sBIC[i] <- which.max(res.fa$sBIC)
  res$model_BIC[i] <- which.max(res.fa$BIC)

#WBIC部分は非掲載とした
}

summary(as.factor(res$model_sBIC))
summary(as.factor(res$model_BIC))

WBIC含めた結果は以下の通りとなった

手法＼推定された因子数 1 2 3 (true) 4 5
sBIC 0 0 80 19 1
WBIC 0 91 8 1 0
BIC 0 0 99 1 0

ということでWBICよりsBICの性能が高いことが分かる。
（WBICの結果の妥当性はやや不安が残るが、sBIC論文でのReduced Rank RegressionでのsBIC vs WBIC や他の因子分析のWBICの結果の報告とも傾向があっているので、傾向としては間違っていないと思われる）

しかし結果を見てみると、これら2つの特異モデル用の情報量規準よりBICの性能が非常に高い！

これについてはJRSS-BのsBIC論文のDiscussionでRafteryがLopes and West (2004)を引用して言及している通り "Good performance of standard BIC was also reported for factor analysis" なので、因子分析については実務的にはBICで良いのかもしれない。

sBICがそもそも何ものなのか、何をついて計算しているものなのかについてはまた今度書こうと思います。

参考文献
Drton M and Plummer M., A Bayesian information criterion for singular models, Journal of the Royal Statistical Society: Series B 2017, 79, 323--380