以下に述べる2重の意味において、今回の状況では「(-2)*残差対数疑似尤度」を比較することはできません。
「一般化線形混合モデル」プラットフォームにおいて、応答の分布として負の二項分布を指定した場合には、「(-2)*残差対数疑似尤度」は、元データに対してではなく、変数変換した疑似応答に対して、対応する重みで重み付けて(正規分布を仮定した線形モデルを)推定したときの、残差対数尤度になっています。
(1) まず、一般的に、混合モデル(正規分布を仮定した場合の線形混合モデル)において、固定効果が異なるモデルの残差対数尤度を比較することはできません。残差対数尤度は、モデルの残差空間だけに限定した尤度になっています。モデルの固定効果を変更したら残差空間も変化するため(喩えるならば、モデルを大きくすると実質的な標本サイズが小さくなるため)、残差対数尤度を比較できません。
(このため、「混合モデル」プラットフォームでは、残差対数尤度だけではなく、「(-2)*対数尤度」(対数尤度の-2倍)も出力しています。)
(2) 今回の場合の残差対数尤度は、元の応答値に対する残差対数尤度ではなくて、それらを変換した疑似応答、および、対応する重みをもつデータに対する(正規分布を仮定した)残差対数尤度になっています。一般的に、異なる変数変換したデータに対する尤度を比較することはできません。たとえば、元データに対する正規分布の尤度と、対数データに対する正規分布の尤度をナイーブに比較することはできません(変数変換した影響を考慮する必要があります)。
もし、変量効果が1つもなくて、固定効果だけの負の二項分布をあてはめたい場合には、「一般化回帰」のほうを用いたほうがいいと思うのですが、負の二項分布の確率分布を"真面目"に扱う場合には、応答変数のデータは整数値でなければいけません。「一般化回帰」プラットフォームにおいても、応答変数のデータは0以上の整数値である必要があります。
直観的な考えですが、対数尤度らしきものを計算するには、以下のような考えがあると思います。
(以下はあくまで私個人による直観的な考えであり、広く使われているわけではありません。)
* もし、平均をとっているデータの個数がすべて同じならば、「平均」が負の二項分布に従っていると考えるのではなく、「合計」が負の二項分布に従ているとして、モデルを変更する。(そもそも負の二項分布の確率関数は、"真面目に"取り扱ったら0以上の整数値しか取りません。一般化線形モデルの疑似尤度による方法では、平均関数と分散関数さえ指定すればいいので、尤度を"真面目に"指定しないでも大丈夫になっています。)
* 予測値をデータテーブルに出力して、自分で何かしらの"対数尤度"らしきものを計算する。
* Wald統計量によって、"対数尤度"らしきものの相対的な大きさを近似する。
今回の場合においてどのように「残差対数尤度」が計算されているかのイメージを伝えるために、切片だけのモデルで計算した例を以下に示します。[ファイル]→[新規作成]→[JSLスクリプト]を選択した後に、呼び出されたウィンドウにおいて以下のコードをコピペし、右クリックして[スクリプトの実行]を選択してください。
// 切片なしのモデルを例としたおもちゃな例
Names Default To Here(1);
nsim = 1000;
y = 5*Negative Binomial Quantile(5, 4, ((1::nsim)`-0.5)/nsim);
// 切片しかない場合には、反復の必要がなく、問題がかなり簡単になる
s = Stddev(y);
my = Mean(y);
alpha = Max(0, (s^2-my)/my^2);
w = my:/(1+alpha * my);
py = Log(my) + (y - my)/my;
Show(alpha, Log(my));
dt = As Table(y||py||J(NRow(y),1,w), <<Column Names({"y", "疑似y", "重み"}));
New Window("Test",
H List Box(
dt << Fit Model(
Y( :y ),
Effects,
Personality( "Generalized Linear Mixed Model" ),
Convergence Limit( 1e-12 ),
Generalized Distribution( "Negative Binomial" ),
Run( Fit( Random Effects Covariance Parameter Estimates( 0 ) ) )
);
,
dt << Fit Model(
Weight( :重み ),
Y( :疑似y ),
Effects,
Personality( "Mixed Model" ),
Run( Random Effects Covariance Parameter Estimates( 0 ) )
);
));
左側が「一般化線形混合モデル」の結果です。右側は、元の応答値を変換した「疑似Y」に対して、それに対応した重みを使って、正規線形混合モデルを推定した結果です。両者の残差尤度関数が一致しているのを確認できると思います。
[注意書き]
上記の文章を作成するにあたり、生成AIで確認した情報を利用しました。次の点を確認・調査するのに生成AIを用いました。
[1] 負の二項分布をIRWLS(反復重み付き法)によってあてはめるためのアルゴリズムの確認
[2] 負の二項分布に従うデータを定数場合したら、負の二項分布には属さないようになることの確認
[3] 疑似尤度法によるモデル推定での変数選択の方法についての文献調査の依頼
Yusuke Ono (Senior Tester at JMP Japan)