要約：線形混合効果モデルでは標準的なAICが根拠を失うので、条件付きAIC(cAIC、AICcとは別物なので注意)が提案されている。Rのlme4パッケージのためにcAICを計算する関数を作製した。
Summary: While Linear Mixed Model have been used in various area, standard Akaike Infermation Criteria might lost its theoritical justification. Here I coded conditional Akaike Information Criteria, proposed by Hodges et al.
線形混合効果モデル（LMM）は素朴な階層ベイズモデルの一種で、生態学、疫学、社会学など、観測データに基づく統計分析を要する分野で広く普及している。しかし、ソフトウェアで簡単に計算できる割に、実はLMMのパラメータ推定法やモデルの評価基準が複雑である。パラメータやモデルの評価は、研究結果を左右するだけに慎重に扱う必要がある。
例えばLMMでは標準的なp値の算出法に根拠がない。例えばRのlme4パッケージでは、作者の配慮によって推定結果にp値を表記されないようになっている。（ただしネットで探すと、p値を算出できるようにするための改造プログラムが出回っている）。
一方で、標準的なAICに関して簡単に計算できる仕様なので多くの混乱を呼んでいる。
例えば、
a<-lmer(y~x1+x2+(1+x1|x3),data=data)
のようにモデルを作り
aic(a)
とコマンドを打つだけで問題なくAICが計算される。しかし、実はLMMで本当にAICを使って良いのか、というと結構論争がある。
##以下、しばらくテクニカルな話題なので結果だけ知りたい人はしばらく読み飛ばしてほしい##
まず混合効果モデルの場合、自由パラメータの数を”どう数えるか”が問題となる。ランダム効果を司るパラメータは、固定効果の切片や偏回帰係数と異なり、正規分布に従って確率的に揺らぐパラメータの分散を定めるものだから、普通のパラメータとは身分が異なる（いわゆるハイパーパラメータ）。
ところでAICにおける罰則項を根拠付けていたのは「期待対数尤度を最大対数尤度で推測した時に生じるバイアスを解析的に解いていくと、その期待値が（一定の前提のもとで）推定されるパラメータの数に一致する」というものであった。
(https://ameblo.jp/yusaku-ohkubo/entry-12257909424.html)
この過程では最尤推定量の漸近正規性などを利用しているが、ハイパーパラメータについても成立するとは限らない。従って、固定効果のパラメータと変量効果のパラメータを同じ基準で数えてしまうとまずい。
次に、最大対数尤度をどう定義するか、という問題がある。つまり、ハイパーパラメータを推定するときの不確実性について、事後分布で重み付平均をした尤度を使うのか、ハイパーパラメータは尤度を最大化する値に固定してしまうのか、という問題である。これは、最下層レベルでの予測を最適化したいのか、階層の一つ上のレベルで最適化したいかによっても変わってくる。以下では、ハイパーパラメータを周辺尤度を最大化する値に固定する場合に有効な方法を紹介する。
##テクニカルな話題終わり##
Hodges and Sargent (2001)は、線形混合効果モデルのために一般化されたcAICを提案した。
これは、通常のAICに出てくる自由パラメータ数をハット行列の対角成分の和（つまりトレース）に置き換えたものだ。ただしハット行列とは、目的変数のベクトルyから、説明変数と最尤推定値で条件づけられたyの期待値への写像を成す行列だ。（なにやら難しそうに聞こえるが、要は「実際のデータプロットと、LMMで作られた回帰線の”ズレ”」をn個のサンプル全てについて計算し格納した行列だ。）
cAICでは、おなじみの「自由パラメータの数」という要素が見当たらないので困惑する人もいるかもしれない。しかし、そもそもAICを「モデルの複雑さに対する罰則」
以下で、Rのソースコードを記載する。ただしLMMオブジェクトの作成時にREML=FALSEを指定しないと計算できないようになっている。
cAIC<-function(object){
H<-hatvalues(object,fullHatmatrix=TRUE) ##Hat行列の準備
tmp<-sum(diag(H)) ##Hat行列のトレースを計算
return(deviance(object)+2*tmp) ##-2loglik+2*tr(H)
}
##使用例
a<-lme4(y~x1+x2+(1+x1|x3),data=data, REML=FALSE)
cAIC(a)
ナイーブに計算されたAICと比較してみよう。混合効果の影響量によっては、ずいぶん異なる値が出てくるはずだ。