またもやマルチレベルモデルの話。
最近統計の話ばっかりだなぁ。まぁいいや。
今回はマルチレベル（特に階層線形モデリング）を実行するときの、設定の注意点についてです。
さて、HLM6やHLM7を使うと、fixed effectの結果のところに、ロバスト標準誤差というのが出力されます。
このロバスト標準誤差ってなんぞや？って思ったことあると思うのですがそういう話です。
あと、HLM7、SPSS、SAS、Mplus、そしてRなどで階層線形モデリングはできますが、それらのソフトウェアの設定の仕方で結果がいろいろ変ります。そのことについても触れておきます。
最後に、ベイズ推定のことについても書いています。ベイズ、オススメです。
大半の人には興味ない話だと思うのですが、興味ある人は続きをどうぞ。
◆ロバスト標準誤差ってなんぞや？
HLMに限らず、誤差やランダム効果の分布には正規分布を仮定しているわけです。
もしこのモデルの仮定が正しければ、ロバストじゃないほうを参照しても問題ないのですが、仮にこれらの正規性が満たされてない場合、つまりモデルの仮定がうまくフィットしてない場合は、パラメータの標準誤差は歪んで推定されます。
ロバスト標準誤差は、こういったモデルの仮定にフィットしてない場合にも適切な検定が行えるように推定がされています。つまり、モデルの仮定の違反に対して頑厳な（そのままですが）わけです。
一般に、ロバスト標準誤差は、普通に最尤法で推定される標準誤差よりも大きく推定されます。なのでモデルにフィットしてなくてもTypeⅠエラーを起こしにくいようになってるんだと思います。（訂正：常にそうとは限りませんでした。）ただ、このロバスト標準誤差は妥当な推定にはある程度のサンプルサイズを必要とするらしく、100以上あるといいようです。逆に言えば、それ以下の場合は怪しい、と。
結果として、どっちを参照したらいいのか？っていう話なんですが、HLMを使うようなデータはサンプルサイズが大きい場合が多いので、ロバストな方を見とけばいいんじゃないかと思います。一応、データの正規性などをよく確認しておく必要はあるとは思いますが。
◆ソフトウェア間で検定結果を一致させるためには？
最近HLMはいろんなソフトウェアで実行できるようになりました。しかし、推定方法や標準誤差の計算の違いによって、p値はいろいろ変わります。ソフトウェア間で結果が変わると困るので、とりあえず一致させる方法をまとめておきます。ただ、完全に一致はしません。今回挙げるのは、HLM7、SPSS、SAS、Mplus、Rです。
まず気を付けるべきは、以下の点です。
- 欠損値の処理方法： リストワイズ削除か否か
- パラメータの推定方法： MLかREMLか
- ロバスト標準誤差か、否か
- 自由度の補正方法： Satterthwaitの補正の有無
まず1の欠損値の方法は、意外と盲点です。
HLMはMDMファイルを作るときに欠損値を省くか、分析を実行するときに省くかを選べます。モデルに投入する変数によって後者は使うサンプル数が変わるので、注意が必要です。
また、MplusなどのSEM用ソフトウェアの場合、独立変数の分散を推定すると、完全情報最尤法によって推定するサンプルが変わったりします。HLMと一致させたいときは、独立変数の分散や 共分散は推定しないようにすることをお勧めします。
サンプルサイズは、基本従属変数のサンプルサイズで決まります。なので、それが一致していなければ当然結果は変わります。
SASとSPSS,Rは汎用ソフトウェアなので、それぞれ欠損値の扱いを同じにすれば結果は同じです。
２の推定方法について
HLMやSPSS,SASそしてRのlme関数で実行すると、デフォルトは制限付き最尤法（REML）という方法で推定します。これは、分散成分の推定が普通の最尤法よりも妥当といわれています（たぶん不偏推定量になる）。大サンプルだと結果はおそらく一致しますが、レベル2のサンプルが小さい場合などはREMLのほうがいいでしょう。一般に、MLのほうがREMLよりも分散パラメータが小さく推定されます。
どのソフトウェアでも最尤法にすることができます。なお、MplusはREMLは選択できません。なのでMplusはほかのソフトウェアと結果を一致させることができないことがあります。
３．標準誤差の推定方法
3は上で触れたロバスト標準誤差の推定。HLMはデフォルトで出力、SASはEmpiricalオプションを加えれば出力、SPSSは出力しません（ただし、一般化混合モデルを使えば一応出せる）。Mplusの場合はMLRで推定すればOKです。Rのlme関数は、そういったオプションはないようですが、自力でなんとかすることはできるかもしれません。標準誤差をロバストにしても、点推定値自体は変わりません。
４．自由度の補正
４の自由度の補正ですが、これが一番厄介です。なぜなら、ソフトウェア間でうまく一致させることができないからです。HLMはレベル2はレベル2のサンプルサイズ、レベル1は全体のサンプルからレベル2のサンプルを引いたものを自由度に使います（また、固定効果の推定分だけ自由度は減る）。ただ、ランダム係数を推定するとレベル1のパラメータもレベル2と同じサンプルサイズになるようです。
Mplusは固定効果を増やしても自由度は減りません。基本的にSEMは検定方法にZ値を使うので、t値を使うHLMやSAS、SPSSとはそもそも方法が違うわけです。SPSSは自動的にSatterthwaitの補正を行います。選択の余地はありません。SASはいろんな補正方法を選べます。HLMと一致させたい場合は、ddfmオプションでbwと指定します（厳密には一致しないがほぼ同じ）。
なぜHLM7などはt検定なのか？という疑問があると思いますが、おそらく制限付き最尤法を使うのがデフォルトだからではないかと思います。
以上、ソフトウェア間でいろいろ設定が違うので、完全に結果を一致させることは困難です。SASとSPSSは同じプロシージャを使っているので完全に一致させることができますが、MplusとHLMは結構結果が異なります。どれがいいんだ、っていう話は一概には言えません。ただ、階層線形モデリングをやるなら今のところHLM7が一番選択肢が多くていいんじゃないかと思います。
◆ベイズの勧め
いろいろ推定法とか標準誤差の計算方法あるみたいだけど、何がベストなの？っていう疑問があるかと思います。
階層線形モデリングみたいに複雑なモデルを計算する場合、従来の最尤法だけではいろいろ問題が出てくるので、それをカバーするためのオプションがたくさん考案されてきたわけです。HLM7はこれらにすべて対応しているので、Bestに近いBetterだと思います。
ただ、そもそも階層線形モデリングという方法にも限界があって、マルチレベルSEMなどの方法に頼りたくなる場合もあるわけです。そうすると、制限付き最尤法が使えなかったり、まぁいろいろ制約がでてきます。
そこで個人的にオススメしたいのはベイズ推定です。最尤法は漸近理論に基づいてるので、サンプルサイズが大きくなるほどよい推定をする、という望ましい性質は持ってはいるものの、小サンプルに対しては怪しい。一方、ベイズ推定は事前分布をうまく活用できれば、小サンプルサイズでも妥当な推定ができます。それに、パラメータの正規性を仮定しないので、最尤法がうまく標準誤差を推定できないようなデータでも、信頼区間を用いて適切に有意性検定を行うことができます。
それ以外にも、マルチレベルモデルとベイズ推定の相性の良さはあります。一般に、レベル2（Betweenレベル）の分散の推定値は正規分布しません。それを最尤法で無理やり推定すると、Wald検定結果は非常に怪しいです。正規分布しにくい理由はいろいろあると思いますが、一つに分散が定義上、負の数を含まないパラメータだから、というのがあると思います。それにもかかわらず最尤法は余裕で0より小さい（つまり負の）分散を推定します。マルチレベルSEMで不適解がよく出るのは、これが原因の一つになっています。
ベイズ推定のMCMCを使うと、パラメータが正規分布してないような場合も無理に標準誤差を使って検定をする必要はありません。それに、非負の制約を簡単にかけられるので、データが少ない場合に引き起こされる類の不適解を回避することができます。これは小さいサンプルサイズを対象としたマルチレベルモデルを分析するときに、とってもありがたい性質です。
追記：ベイズ推定についてのWSで，MCMCでマルチレベルモデルという内容で発表しました。
そういうわけでベイズは個人的にいいなぁって思うんですが、いかんせんマルチレベルだけでもマニアック扱いされるのに、そこにさらにベイズとか言ってたら、オーディエンスがついてこないので、なんとかならんもんかと思う次第です。
そんなわけで、マルチレベルモデルのマニアックな話でした。