線形モデルと線形混合モデルにおけるAIC

はじめに

千葉大学・株式会社Nospareの川久保です．
今回は線形モデルと線形混合モデルにおけるAIC（赤池情報量規準，Akaike information criterion）型の情報量規準について紹介します．この話題については約50年にわたる膨大な研究の蓄積があり，解説文も日本語英語問わず数多く存在しますが，私が統計学の勉強を始めて最初に感銘を受けたテーマなので，頑張って書いていきたいと思います．

AICの一般論

AICはそれぞれの候補モデルについて計算され，小さい値をとるモデルを良いモデルと判断します．他の情報量規準についても同様で，小さい値が良いモデルです．

設定

まずは特定のモデルのクラスを想定せずに，一般論としてAICの考え方について説明します．観測データ$y = (y_1,\dots,y_n)^\top$の真の密度関数を$g(y)$とします．$g(\cdot)$をデータ解析者は知らないので，$f(y \mid \theta)$という密度をもつパラメトリックモデルを当てはめることを考えます．AICは$y$の独立複製$\tilde{y}$の予測パフォーマンスが良くなるようなモデルを選ぶことを目指しています．$\tilde{y}$が$y$の独立複製であるとは，$y$と$\tilde{y}$が独立に同じ$g(\cdot)$から生成されることを意味します．$y$と同じ分布にしたがう将来の値$\tilde{y}$を予測する，といったイメージです．

KLリスク

さて，$\theta$の（最尤）推定量を$\hat{\theta}$とします．最尤推定量は観測データ$y$の関数と見ることができるので，明示的に $\hat{\theta}(y)$ と書いておきます．$\hat{\theta}(y)$を用いた$\tilde{y}$の予測密度として，$f(\tilde{y} \mid \hat{\theta}(y))$ が考えられます．これは，モデルとしての密度 $f(\tilde{y} \mid \theta)$ において，パラメータ$\theta$に最尤推定量（値）をプラグインしたものなので，プラグイン予測密度と呼びます．このプラグイン予測密度（モデルからデータ分析者が作ったもの）が，$\tilde{y}$の真の密度$g(\tilde{y})$にどれだけ近いかによって，モデルの良さを測ります．その近さとして，以下のKullback–Leibler（カルバック・ライブラー，KL）ダイバージェンスを考えます:
$$
\int \log \Big\{ {g(\tilde{y}) \over f(\tilde{y} \mid \hat{\theta}(y))} \Big\} g(\tilde{y})d\tilde{y}
$$
この量は，最尤推定量がデータ$y$の関数であることから，データ$y$についてのランダムネスを持っています．つまり，観測データ$y$がどのような値が実現するか（標本として得られるか）によって，値が異なってきます．そこでデータ$y$の真の分布（密度）$g(y)$で期待値をとります:
$$
\underbrace{\int \bigg[ \overbrace{ \int \log \Big\{ {g(\tilde{y}) \over f(\tilde{y} \mid \hat{\theta}(y))} \Big\} g(\tilde{y})d\tilde{y} }^{損失} \bigg] g(y)dy }_{リスク}
$$
この量をKLリスクと呼ぶこととしますが，それは統計的推測決定理論における損失関数とリスク関数の考え方に通じるからです．上式のKLリスクのうち，モデルに関係する項のみを取り出し2倍した量をAkaike information (AI) とします:
$$
\mathrm{AI} = -2 \iint \log \{ f(\tilde{y} \mid \hat{\theta}(y)) \} g(\tilde{y}) g(y)d\tilde{y} dy
$$
このAIの（漸近）不偏推定量がAICで，よく知られている以下の形になります．
$$
\mathrm{AIC} = -2 \log \{ f(y \mid \hat{\theta}(y)) \} + 2p
$$
第1項は最大対数尤度（対数尤度関数を最尤推定値で評価した量）の-2倍，第2項の$p$はモデルの未知パラメータ数（$\theta$の次元）です．AIC値が小さいほど良いモデルと評価されますが，それはAICがリスクの推定量だからです．

第2項の解釈

AICの第2項はバイアス補正項またはペナルティ項と呼ばれますが，2つの名称は由来が違いますので，それぞれについて簡単に解説します．

バイアス補正項

AIを最大対数尤度の-2倍（すなわちAICの第1項）で推定すると，下方バイアス，すなわちAIというリスクを過小に推定してしまいます．直感的には，最大対数尤度というのはデータを見てから，それに最も当てはまる（尤もらしい）パラメータの値を代入しているので（後出しジャンケンとよく例えられます），データに当てはめ過ぎているわけです．そのバイアスを評価すると漸近的に$-2p$になることが知られていて，バイアス補正項として$2p$を加えてあげると，$n\to\infty$ でAICの期待値がAIに等しくなります（AICはAIの漸近不偏推定量）．

ペナルティ項

上でも書きましたが，AICの第1項はモデルの当てはまり（の悪さ）を表していますが，モデルを複雑にすれば（パラメータの数を増やせば），第1項を小さくしていくことができます．ところがパラメータの数が多い複雑なモデルは推定の効率が落ち，一般に予測精度も悪くなります．パラメータの数（$p$）が大きくなれば，AICの第2項の$2p$は当然大きくなるので，この項をモデルの複雑さに対するペナルティ項と解釈できるわけです．

BICとの比較

AICとならぶ代表的な情報量規準として，BIC（ベイジアン情報量規準）があります:
$$
\mathrm{BIC} = -2 \log \{ f(y \mid \hat{\theta}(y)) \} + \log(n)\cdot p
$$
第1項はAICと全く同じです．第2項のペナルティについては，$n \geq 8$ のとき $2 < \log(n)$ なので（$n$はサンプルサイズ），$n \geq 8$ においては（ほとんどの実データではそうだと思いますが）BICの方がAICよりも強くペナルティがかかった規準と言えます．つまりBICの方がAICよりも小さいモデル（パラメータの数が少ないモデル）を選びやすいです．

線形モデルの例

以下の線形モデルを考えます:
$$
y = X\beta + \varepsilon, \quad \varepsilon \sim \mathrm{N}(0,\sigma^2 I_n)
$$
説明変数の数（$\beta$の次元）を$k$とすると，このモデルにおける未知パラメータは$\beta$と$\sigma^2$なので，$p = k+1$です．最大対数尤度を整理すると，以下が得られます:
$$
\mathrm{AIC} = n \{ \log(2\pi\hat{\sigma}^2) + 1 \} + 2(k+1),
$$
ただし，$\hat{\sigma}^2 = n^{-1}(y-X\hat{\beta})^\top(y-X\hat{\beta}),\ \hat{\beta} = (X^\top X)^{-1}X^\top y$ です．
Rのlm()関数を使って線形モデルを推定した場合，以下のようにAIC()関数を使って簡単にAIC値を計算できます．

res <- lm(y ~ x)
AIC(res)

線形混合モデルにおけるAIC

線形モデルの拡張として，以下のような線形混合モデルを考えます:

y = X\beta + Zb + \varepsilon, \quad b \sim \mathrm{N}(0,G), \quad \varepsilon \sim \mathrm{N}(0,R)

$X\beta$の項を固定効果（fixed effect），$Zb$の項を変量効果（random effect）と呼び，2つの効果が"混ざった"モデルということで，線形混合モデルと呼ばれています．$Z$のデザインによって，個人効果や空間効果などを自由にモデリングできます（線形混合モデルの応用事例などは別にいつか記事を書きたいと思っています）．以後の説明のために，いくつか記法を整理しておきます．

\begin{align}
f(y \mid b,\theta) &\sim \mathrm{N}(X\beta + Zb, R), \\
\pi(b \mid \theta) &\sim \mathrm{N}(0,G), \\
f_\pi(y \mid \theta) &= \int f(y \mid b)\pi(b)db \sim \mathrm{N}(X\beta, \Sigma), \quad \Sigma = ZGZ^\top + R
\end{align}

$\theta$はこのモデルにおける未知パラメータで，共分散行列$G$と$R$に含まれる未知パラメータと，回帰係数$\beta$の集まりです．

marginal AIC

このモデルに対するAICは，変量効果の取り扱い方によって2通り考えられます．まず1つが，変量効果$b$を積分消去した$y$の周辺密度（marginal density）$f_\pi(y \mid \theta)$に対して，AICを定義する方法です．これをmarginal AICと呼びます:
$$
\mathrm{mAIC} = -2\log \{ f_\pi(y \mid \hat{\theta}) \} + 2p,
$$
やはりここでも$p$は未知パラメータの数，すなわち$\theta$の次元です．
mAICは以下のKLリスクの推定を目指しています:
$$
\underbrace{ \int \bigg[ \overbrace{ \int \log \Big\{ {g_\pi(\tilde{y}) \over f_\pi(\tilde{y} \mid \hat{\theta})} \Big\} g_\pi(\tilde{y}) d\tilde{y} }^{損失} \bigg] g_\pi(y)dy }_{リスク}
$$

ただし$g_\pi(\cdot)$は，$y$ないし$\tilde{y}$の真の周辺密度です．

conditional AIC

Vaida and Blanchard (2005) は，変量効果$b$を所与とした$y$の条件付密度（conditional density）$f(y \mid b,\theta)$をもとに構成した，conditional AIC (cAIC) を提案しました:
$$
\mathrm{cAIC} = -2\log \{ f(y \mid \hat{b},\hat{\theta}) \} + 2\rho
$$
第2項（バイアス補正項，ペナルティ項）における$\rho$（ロー，$p$ではないですよ！）は，線形混合モデルにおけるモデルの複雑さとして提案されていたeffective degrees of freedomと呼ばれるものです．cAICは以下のKLリスクの推定を目指しています:
$$
\underbrace{ \iint \bigg[ \overbrace{ \int \log \Big\{ {g(\tilde{y} \mid b) \over f(\tilde{y} \mid \hat{b},\hat{\theta})} \Big\} g(\tilde{y} \mid b) dy }^{損失} \bigg] g(y \mid b) \pi(b) dy db }_{リスク}
$$
ただし$g(\cdot \mid b)$は，$y$ないし$\tilde{y}$の$b$を所与とした真の条件付密度です．

mAICとcAICの違い

注意が必要なのは，mAICにおける$\tilde{y}$と，cAICにおける$\tilde{y}$は意味が違うということです．mAICにおける観測モデル（$y$に対するモデル）と予測モデル（$\tilde{y}$に対するモデル）はそれぞれ以下のように表されます:

\begin{align}
観測モデル: & \quad y = X\beta + u, \quad u \sim \mathrm{N}(0,\Sigma), \\
予測モデル: & \quad \tilde{y} = X\beta + \tilde{u}, \quad \tilde{u} \sim \mathrm{N}(0,\Sigma)
\end{align}

つまり線形混合モデルにおける変量効果$Zb$と誤差$\varepsilon$をまとめた$u = Zb + \varepsilon$に対して，その独立複製$\tilde{u}$を予測モデルで想定しています．一方cAICにおける観測モデルと予測モデルは以下です:

\begin{align}
観測モデル: & \quad y = X\beta + Zb + \varepsilon, \quad \varepsilon \sim \mathrm{N}(0,R) \\
予測モデル: & \quad \tilde{y} = X\beta + Zb + \tilde{\varepsilon}, \quad \tilde{\varepsilon} \sim \mathrm{N}(0,R)
\end{align}

誤差$\varepsilon$についてのみ，その独立複製$\tilde{\varepsilon}$を予測モデルでは考えており，変量効果$Zb$は観測モデルと予測モデルで共通です．線形混合モデルの応用では，変量効果は以下のような効果としてモデリングすることが考えられます．

• 個人効果（経時データ，パネルデータ）
• 地域効果（空間データ）

つまりデータに何らかのクラスターを想定し，そのクラスターの効果を変量効果でモデリングしています．こうしたモデリングを行ったときは，同じクラスター内での将来の値を予測したい，といった要請が多く，その際はcAICが想定している「観測モデルと予測モデルの変量効果は同じ」という設定の方が自然です．

RのcAIC4パッケージ

RにはcAIC4というパッケージでcAICが実装されています．パッケージ内のcAIC()関数によりcAIC値を計算できますが，その引数にはlme4というパッケージ内の関数で返されたオブジェクトが必要になります．以下，cAIC()関数のヘルプ内のExampleを引用します．

library(cAIC4)
b <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
cAIC(b)

lmer()はlme4パッケージにおいて線形混合モデルを当てはめる関数，sleepstudyはlme4内のデータセットで，睡眠不足の日数Daysが反応時間Reactionに与える影響を調べたデータセットです（結構とんでもない実験をやっているんですね…！）．

まとめ

今回はAICの一般論と，線形混合モデルに対するcAICについて説明しました．AICを含む情報量規準は以下の本に詳しいです（私もこの本で最初勉強しました）:

小西貞則・北川源四郎「情報量規準」

cAICが初めて提案されたのは以下の論文です:

• Vaida, F. and Blanchard, S. (2005). Conditional Akaike information for mixed-effects models. Biometrika, 92, 351–370. DOI: 10.1093/biomet/92.2.351

また，私たちはcAICを改良した情報量規準をいくつか論文で発表しています．

• Kawakubo, Y. and Kubokawa, T. (2014). Modified conditional AIC in linear mixed models. Journal of Multivariate Analysis, 129, 44–56. DOI: 10.1016/j.jmva.2014.03.017

• Kawakubo, Y., Sugasawa, S. and Kubokawa, T. (2018). Conditional Akaike information under covariate shift with application to small area estimation. Canadian Journal of Statistics, 46, 316–335. DOI: 10.1002/cjs.11354, arXiv

こちらの話題についても，今後記事として紹介できればと考えています．

株式会社Nospareには，統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospare までお問い合わせください．