StanなMCMCで探索的因子分析やってみた

この記事はStan Advent Calender 2016にエントリーしたものです。

やっていることは，ここ数年関わっている探索的因子分析をベイズ推定するっていう話なのですが，最初に言っておきますと，探索的因子分析ってベイズモデリングに向いてないと思うのですよw

ベイズモデリングが何か，という話は一言では終わらないですが，私は「データ生成メカニズムをデザインするモデリング」と理解しています。なので，認知心理学的なモデルなんかにはもってこいだと。そもそも心理学は個人差をどう考えるか，というような世界でもあるので，階層的なモデリングのためにとてもいいと思うのですね。

一方，探索的因子分析はその名の通り「探索的に」データを捉えるわけですから，データがどういう構造しているかわからない，というのはそもそもちょっと違う。もちろん，n因子構造だと仮定してモデルを書くと，ということになるんですが，CFAならともかくEFAだとまあありとあらゆる係数を推定していくことになるのです。なので，ちょっと大変だなあ，と書き終わってから思うたわけです。

 

さてさて，本題です。
因子分析というのは，$$N$$人の個人$$i$$が$$M$$個ある項目$$j$$に反応したデータセット$$x_{ij}$$(ただし$$N \ll M$$,つまり項目数より個人の数が十分に多い)について，$$k$$個の因子という潜在変数を仮定して，因子と項目との関係を表す因子負荷量$$a_{jk}$$と，因子と個人との関係を表す因子得点$$f_{ik}$$との積和の形で表すモデルです。

$$ x_{ij} = \mu_j + \sum_{k=1} a_{jk}f_{ik} + e_{ij}$$

ここの$$e_{ij}$$は誤差項ね。これをベイズ的に考えるときは，基本的に潜在変数が正規分布に従うことを仮定する，というだけなので，それをそのまま素直にモデリングしてあげれば一丁あがりです。

ただ，探索的にやる場合は因子負荷行列に制約をかけてあげる必要があります。因子負荷行列は一般に，直交行列になるような制約を置いてあげるものですが，SEM的には$$a_{jk}=0$$(ただし j < k の場合)という制約をかけることになります(豊田先生の「共分散構造分析【技術編】」を参考）。これをコード化したら次のようになります。

data{
  int<lower=0> A; // number of factors
  int<lower=0> M; // number of items
  int<lower=0> N; // number of observations
  int<lower=0> D;
  vector[M] r[N];
}

parameters {
  vector[A] F[N];                 //factor scores
  real<lower=0> e[M];       //SD of unique factor
  vector[D] ELT;                  //Elements of the low triangle matrix(factor loadings)
}

transformed parameters { 
  matrix[M,A] L; //constrained factor loading matrix
  {
  int index;
  index = 1;      
  for(m in 1:M)
    for(a in 1:A)
      if(a > m){
        L[m,a] = 0;
      }else{
        L[m,a] = ELT[index];
        index = index + 1;
      }
  }
}

model{
  ELT ~ normal(0.1,100);
  e ~ cauchy(0,5);
  for(n in 1:N)
    F[N] ~ normal(0,1);
  
  for(n in 1:N)
    r[n] ~ normal(L*F[n],e);
}

ここのtransformed parametersブロックで面倒なことをしています。整数でindexというのを宣言して置いて，このループの中だけでカウントするインデックスを用意しました。もう少しうまいやり方もあると思うんだけど，とりあえずこれでできます。

 

ただし！このやり方はすごく時間がかかります。時間がかかる理由は，データサイズに比例して推定しなければならないパラメタが多いということ。N人A因子あればN*A分の因子得点を推定しますからね！こりゃめんどうです。

実際，一般的な因子分析は，基本モデルにいくつかの条件（因子得点は標準化得点，誤差の平均はゼロ，誤差因子と共通因子の相関はゼロ)を加えて，

$$ \Sigma = A \Phi A’ + \Delta $$

特に直交因子モデルの場合は$$\Phi$$は単位行列なので

$$ \Sigma = AA’ + \Delta $$

とした分散共分散行列で考えることになります。

Stanで多変量正規分布を推定するのはなぜか大変重たいので，ウィシャート分布から推定することにし，因子負荷行列もこれスキー分解したベクトルで与えてやると，Stanコードは次のように大変短くなります。

data{
  int<lower=0> A; // number of factors
  int<lower=0> M; // number of questions(items)
  int<lower=1> N; // number of subjects
  matrix[N,M] y; // centered dataset
}

transformed data{
  vector[M] mu; // zero means
  matrix[M,M] spd; //sum of products of deviation
  mu = rep_vector(0,M);
  spd = y'*y;
}

parameters {
  cholesky_factor_cov[M,A] L;   //constrained factor loadings
  vector<lower=0>[M] u2;   //variances of unique factor
}

transformed parameters {
  cov_matrix[M] Sig; //estimated covariance matrix
  Sig = tcrossprod(L)+ diag_matrix(u2);
}

model{
  to_vector(L) ~ normal(0,100);
  u2 ~ cauchy(0,5);
  spd ~ wishart(N-1,Sig);
}

データから平均値を取り去った(中心化した)データセットとして与え，それの積和行列がtransformed dataブロックでspdとして作られており，それがwishart分布に従うようになっています。

これでStanはちゃんと推定してくれます。が！これでお話は終わりではありません。実はここから得られるMCMCサンプルの推定値を考えても，0ばっかり！みたいなことになるのです。それには三つほど問題が関係しています。

一つはMCMCサンプルごとに出てくる因子の順番が違う問題。これはサンプリング後に因子の順番をソートして整えることで対応できます。

もう一つはMCMCサンプルごとに出てくる因子負荷量の符号が違うこと。これは因子毎の負荷量の総和が正になるように整えるとか，任意のMCMCサンプルをつかってその因子毎の負荷量の総和の符号に合わせる，といった恣意的な基準を作って対応できます。

もう一つは因子の回転についてです。stanの中では回転作業をしていませんので，得られたMCMCサンプルを全部回転することで対応できます。

 

・・・とここまで書いておわかりのように，サンプルを得た後の作業が結構大変なんですね。実際にやってみると，サンプリングよりその後の回転やソートのほうが時間がかかったりします。

その作業をするコードを載せても良いのですが，それはRの話でStanの話じゃなくなってしまいますから，この記事はここまでです。Rの関数として，psychパッケージのようにfa，と書いたら探索的ベイジアン因子分析してくれるものを現在開発中ですので，興味がある方は個別にご連絡を。

 

それにしても，確率モデルをほぼそのまま書くだけでできちゃうStanって，本当に便利ですよね。みなさんもEnjoy Stan!

参考文献；

[amazonjs asin=”4254126646″ locale=”JP” title=”共分散構造分析 技術編―構造方程式モデリング (統計ライブラリー)”]