GLMM の基礎

GLMM タブで使われている統計理論の背景です。操作方法は GLMM のページを参照してください。

階層データと独立性の問題

統計モデルの多くは、観測が互いに独立であることを仮定しています。しかし現実のデータには、観測がグループに属する階層構造を持つものが多くあります。

• 複数の学校に通う生徒のテスト成績（生徒はそれぞれの学校に所属）
• 複数の病院で治療を受けた患者の回復日数（患者はそれぞれの病院に所属）
• 同じ被験者に対する反復測定（測定は被験者に所属）

同じグループに属する観測は、グループ固有の要因（学校の教育方針、病院の設備、被験者の体質）を共有するため、互いに似た値を取る傾向があります。この相関を無視して GLM で分析すると、モデルの誤特定（misspecification）になります。データの生成過程にグループ構造があるのに、モデルがそれを表現していない状態です。

この誤特定は、推定結果を 2 つの形で損ないます。

標準誤差の過小評価: グループ構造を無視した GLM は、観測間の共分散をゼロと仮定したまま標準誤差を計算します。実際には同じグループ内の観測に正の相関があるため、切片やグループ間で変動する説明変数の係数で標準誤差が過小になり、95% 信頼区間は実際には 95% の被覆確率を持ちません。過小の度合いは説明変数の変動がグループ間とグループ内のどちらにあるかで変わります（ICC の節のデザイン効果を参照）。

固定効果の偏り: グループレベルの交絡要因があると、固定効果の推定に欠落変数バイアス（omitted variable bias）が生じます。たとえば学校の教育方針が勉強時間と成績の両方に影響する場合、学校差を無視した回帰では勉強時間の効果を正しく推定できません。ランダム切片はグループ間の平均レベルの違いをモデル化しますが、その推定はランダム切片が説明変数と無相関であることを仮定しています。未測定のグループレベル交絡が説明変数と相関している場合はこの仮定が崩れ、固定効果の推定は偏ったままです。測定された交絡要因は固定効果として投入する必要があります。

混合モデル（mixed model）は、グループ内の相関をランダム効果として明示的にモデル化することで、標準誤差の過小評価を解消します。あわせて、ばらつきをグループ間とグループ内に分けて推定するため、グループの違いに由来するばらつきの大きさを定量的に把握できます。この分解を要約する指標が後述の ICC です。

固定効果とランダム効果

混合モデルの「混合」は、固定効果とランダム効果の両方を含むことを指します。

固定効果（fixed effects） は、どのグループにも共通に働くと仮定する、説明変数と応答の体系的な関係です。係数は未知の定数としてモデル化されます。たとえば「勉強時間が 1 時間増えるとテスト成績が何点上がるか」という効果は固定効果です。

ランダム効果（random effects） は、直接観測できないグループレベルのばらつきをモデル内で表現するための装置で、変量効果とも呼ばれます。各学校の平均成績は教育方針や生徒層によって異なります。このばらつきは確率的なメカニズムで生まれたものではなく、具体的な要因に起因しますが、それらの要因をすべて説明変数として測定・投入するのは現実的ではありません。ランダム効果はこの「観測できないグループ差」を確率分布で近似的に表現します。

グループ $j$ の差を $u_j$ と書くと、$u_j \sim N(0, \sigma_u^2)$ という仮定は「学校差が正規分布から生成されている」という主張ではありません。グループ間のばらつきの大きさ $\sigma_u^2$ を推定し、それに基づいて縮小推定（後述）を行うための、モデリング上の仮定です。

両者の実質的な違いは、パラメータの数え方に現れます。グループをダミー変数にして固定効果として扱うと、グループの数だけ自由なパラメータが増え、各グループの効果はそれ自体が推定対象になります。ランダム効果では、推定する分散成分は $\sigma_u^2$ の 1 つだけで、各グループの $u_j$ は自由なパラメータではなく、この分散を前提として予測される値です。固定効果/ランダム効果の区別は、グループがランダムに抽出されたかどうかではなく、グループ間の差をこのどちらの形で扱うかの選択です（Gelman, 2005）。

ランダム切片モデル

最も基本的な混合モデルがランダム切片モデルです。各グループに固有の切片（ベースライン）を与えます:

$$g(\mu_i) = x_i'\beta + u_{j[i]}, \quad u_j \sim N(0, \sigma_u^2)$$

$g$ はリンク関数、$x_i'\beta$ は固定効果の線形予測子です。$j[i]$ は観測 $i$ が属するグループを表す添字で、$u_{j[i]}$ はそのグループのランダム切片です。$u_j \sim N(0, \sigma_u^2)$ は、グループごとの値 $u_1, \dots, u_J$ がそれぞれ独立に同じ正規分布に従うことを表します。

直感的には、全グループに共通する回帰直線（$x_i'\beta$）を引き、各グループの直線はそこから $u_j$ だけ上下にずれる、という構造です。正規分布 $N(0, \sigma_u^2)$ の仮定は、このずれの大きさを $\sigma_u^2$ という1つのパラメータで要約するための道具です。各 $u_j$ がどんな値をとるかは、データから推測します（ランダム効果の予測と縮小推定を参照）。

Gaussian 分布族の場合、これは線形混合モデル（LMM）です:

$$y_i = x_i'\beta + u_{j[i]} + \varepsilon_i, \quad u_j \sim N(0, \sigma_u^2), \quad \varepsilon_i \sim N(0, \sigma_e^2)$$

$\sigma_u^2$ はグループ間のばらつき、$\sigma_e^2$ はグループ内（個人レベル）のばらつきです。残差分散 $\sigma_e^2$ は全グループ共通と仮定されます。ランダム切片モデルが表すのはグループ間の平均レベルの違いであり、グループ内のばらつきの大きさがグループごとに異なる構造はモデル化しません。

同じグループの観測が互いに似るという性質は、この式では共分散として表れます。同じグループの 2 つの観測は $u_j$ を共有するため、その共分散は $\sigma_u^2$ です。異なるグループの観測どうしの共分散は 0 です。各観測の分散は $\sigma_u^2 + \sigma_e^2$ なので、同一グループ内の 2 観測の相関は $\sigma_u^2 / (\sigma_u^2 + \sigma_e^2)$ となり、この量が後述の ICC です。

リンク関数が恒等でない場合、係数 $\beta$ の解釈は条件付きになります。$\beta$ が表すのは、同じグループの中で説明変数を動かしたときのリンク尺度上の効果です。$g^{-1}$ が非線形なため、グループを均した集団平均の応答に対する効果は、$\sigma_u^2 > 0$ である限り一般にこれと一致しません。たとえば Binomial + logit で係数から得られるオッズ比は、グループをまたいだ平均的なオッズの比ではなく、同一グループ内での比較として読みます。

パラメータ推定

混合モデルのパラメータ推定は、固定効果 $\beta$ と分散成分（$\sigma_u^2$, $\sigma_e^2$）を同時に推定する必要がある点で、通常の回帰モデルの推定より複雑です。

REML（制限付き最尤法）

最尤法（ML）は分散成分を過小推定する傾向があります。$\beta$ の推定に自由度を消費することを尤度関数が考慮しないためです。分母に $n$ を使った標本分散が分散を過小推定するのと同じ現象で、標本分散では分母を $n-1$ にして補正します。

REML はこの偏りを避けるため、$\beta$ の影響を受けない量だけから分散成分を推定します。応答ベクトル $y$ を $n$ 次元空間の点とみると、固定効果 $X\beta$ が動けるのは計画行列 $X$ の列が張る $p$ 次元の部分空間です。$y$ をその直交補空間へ射影した残差は、$X\beta$ の成分がちょうど消えるため、分布が $\beta$ に依存せず分散成分だけで決まります。REML はこの残差の尤度を最大化します。$\beta$ の推定に自由度を消費しないので、分散の推定は実質 $n - p$ 次元のデータに基づき、偏りが補正されます。標本分散の $n - 1$ は、$X$ が切片のみ（$p = 1$）の場合にあたります。

MIDAS は REML の最大化を、探索を 1 次元に落とした profile REML で行います。分散比 $\theta^2 = \sigma_u^2/\sigma_e^2$ を 1 つ与えるごとに、その下で最適な $\beta$ と $\sigma_e^2$ は閉形式（closed form）で決まります。数値的に探索するのは $\log\theta$ だけです。REML 対数尤度は $\log\theta$ について複数の局所的な最大点を持つことがあるので、MIDAS は探索範囲に等間隔に置いた点で REML 対数尤度を評価し、値が最大の点の両隣で挟んだ区間に黄金分割法を使います。この手順を選んだ理由と、この手順でも最大点を見落とす場合は、Laplace 近似と PIRLS の節で外側ループについて説明します。

Laplace 近似と PIRLS

Gaussian + identity ではランダム効果についての積分を解析的に計算できるので、MIDAS は前節の REML で推定します。それ以外の組み合わせ（Binomial、Poisson、Gamma、Gaussian + log）では、ランダム効果について積分した周辺尤度

$$L(\beta, \sigma_u^2, \phi) = \int \prod_i f(y_i \mid \mu_i, \phi) \cdot f_u(u) \, du$$

を閉形式で表せません。$f$ は選択した分布族における各観測の密度関数（Binomial や Poisson などの離散分布族では確率質量関数）、$f_u$ はランダム効果の密度関数、すなわち $N(0, \sigma_u^2)$ の密度です。$\mu_i$ は $\beta$ と $u$ から決まる観測 $i$ の平均で、$\phi$ は分布族の分散パラメータです。MIDAS は Gaussian + log と Gamma では $\phi$ を推定し、Poisson と Binomial では $\phi = 1$ に固定します。

Laplace 近似は、被積分関数の対数を最大点の周りで 2 次式に近似し、被積分関数を正規分布の密度の形に置き換えてこの積分を計算します。グループあたりの観測数が少ないときと、Binomial でイベントが稀なときに、この近似の精度は下がります。

MIDAS は、ランダム効果の分散を相対共分散パラメータ $\theta$ を使って $\sigma_u^2 = \theta^2\phi$ と表し、Laplace 近似した周辺尤度を次の式で計算します。Gaussian では $\phi = \sigma_e^2$ なので $\theta^2 = \sigma_u^2/\sigma_e^2$ です。ただし Gaussian + log の $\sigma_u$ は線形予測子（対数）の尺度の量なので、この $\theta$ は REML の節の $\theta$ と違って応答の単位に依存します。

$$-2 \log L \approx \sum_j \log(1 + \theta^2 C_j) + \frac{1}{\phi}\left(\sum_i d(y_i, \hat\mu_i) + \sum_j \frac{\hat u_j^2}{\theta^2}\right) + \sum_i a(y_i, \phi)$$

式の $d$ と $a$ は、各観測の密度を $-2 \log f(y \mid \mu, \phi) = d(y, \mu)/\phi + a(y, \phi)$ と分けたときの項です。$d(y, \mu)$ は観測の逸脱度への寄与で、$a(y, \phi)$ は $\mu$ を含まない項です。この分け方では、グループ $j$ の被積分関数の対数が、$-\{\sum_{i \in j} d(y_i, \mu_i) + u_j^2/\theta^2\}/(2\phi)$ に $u_j$ を含まない項を足した形になります。そのため被積分関数を最大にする $u_j$ の値 $\hat u_j$ は、$\theta$ と $\beta$ だけで決まります。$\hat\mu_i$ は $u_j = \hat u_j$ としたときの観測 $i$ の平均です。$C_j$ は、グループ $j$ の各観測について $d$ を線形予測子で 2 回微分して $u_j = \hat u_j$ で評価した値の半分を、グループ内で合計したものです。式の第 1 項は、ランダム効果の分布の分散 $\sigma_u^2$ と、被積分関数を置き換えた正規分布の密度の分散 $\phi/(C_j + 1/\theta^2)$ の比の対数を、グループについて合計したものです。

MIDAS は、この近似した周辺尤度を固定効果 $\beta$、$\theta$、$\phi$ のすべてについて最大化します。$\hat\phi$ は、$\theta$ と $\beta$ を与えるごとに上の式を $\phi$ について最小にする値です。Gaussian + log では $\hat\phi = \{\sum_i d(y_i, \hat\mu_i) + \sum_j \hat u_j^2/\theta^2\}/n$ です。Gamma では、MIDAS は対応する方程式を数値的に解いて $\hat\phi$ を求めます。報告する対数尤度・$\phi$ の推定値・$\sigma_u^2 = \theta^2\phi$ は、すべてこの $\hat\phi$ による値です。Gaussian + log の $\hat\phi$ は残差分散 $\sigma_e^2$ の推定値ですが、Gamma の $\hat\phi$ は応答の条件付き分散 $\phi\mu^2$ の係数で、分散そのものではありません。最大化は数値的な反復によるので、収束した場合に得られるのは局所的な最大点で、大域的な最大点とは限りません。反復が収束しなかった場合、MIDAS は推定値に警告を付けるか、推定値を定められずに当てはめを中止します。どちらになるかは GLMM タブの収束の問題の節を参照してください。

最大化の手順は、Bates et al. (2015) が記述した入れ子の 2 ループ構造にもとづきます。外側ループの最適化手法と、内側ループで固定効果を求める手順は MIDAS 独自の選択です。

• 外側ループ: 近似した周辺尤度を、相対共分散パラメータ $\theta$ の対数について、探索範囲に置いた点での評価と黄金分割法で最大化します
• 内側ループ: $\theta$ を固定し、PIRLS（Penalized IRLS）の解から出発して、近似した周辺尤度を $\beta$ について最大化します

外側ループが黄金分割法の前に複数の点で評価するのは、近似した周辺尤度を $\theta$ 以外のパラメータについて最大化した値が、$\theta$ について複数の局所的な最大点を持つことがあるためです。$\theta = 0$ の境界と $\theta > 0$ の点の両方が局所的な最大点になり、$\theta > 0$ の点のほうが尤度の高いデータがあります。黄金分割法は、最大点が 1 つであることを前提に、2 点での値の比較で探索範囲の片側を捨てます。そのため探索範囲の全体に黄金分割法を使うと、$\theta > 0$ の最大点を含む側を捨てて境界解に止まることがあります。MIDAS は $\log\theta$ の探索範囲に等間隔に置いた点で近似した周辺尤度を評価し、値が最大の点の両隣で挟んだ区間に黄金分割法を使います。黄金分割法は値が最大の点から区間を絞るので、得られる尤度はその点での値を下回りません。尤度の最も高い局所的な最大点がこの区間の外にあると、MIDAS はその最大点を見落とします。区間の中に局所的な最大点が複数あるときも、得られるのはそのうちの 1 つです。

PIRLS は GLM の IRLS にランダム効果のペナルティ項を加えたもので、次のペナルティ付き逸脱度を $(\beta, u)$ について最小化します:

$$D_p = \sum_i d(y_i, \mu_i) + \sum_j \frac{u_j^2}{\theta^2}$$

$D_p$ の第 1 項はデータへの適合度、第 2 項は $u_j$ が 0 から離れるほど大きくなるペナルティです。第 2 項は $u_j$ の分布 $N(0, \sigma_u^2)$ の密度の対数に由来し、$\theta$ が小さいほど重くなります。このペナルティにより、データの少ないグループの予測値は全体平均に引き寄せられます（後述の縮小推定を参照）。

PIRLS の解の $\beta$ は、近似した周辺尤度の最大点とは一般に一致しません。PIRLS は $D_p$ だけを最小にしますが、近似した周辺尤度の式の $\sum_j \log(1 + \theta^2 C_j)$ も $C_j$ を通じて $\beta$ に依存するためです。MIDAS は PIRLS の解を出発点として、この項を含めた近似周辺尤度を $\beta$ について最大化します。$\beta$ を動かすたびに、$D_p$ を $u_j$ について最小にする値として $\hat u_j$ を求め直します。

境界解（singular fit）

分散 $\sigma_u^2$ の推定には 0 という下限があります。どちらの推定経路でも、最適化の解がこの下限に張り付いた状態を singular fit と呼びます。境界解は、グループ間の変動が残差の変動に比べて小さいデータで起きます。$\log\theta$ の探索が $\theta = 0$ の境界の局所的な最大点に止まり、境界から離れた位置にある尤度のより高い最大点を見落とした場合にも起きます。探索が最大点を見落とす条件は Laplace 近似と PIRLS の節で説明します。

境界解は、推定結果の読み方を変えます。$\hat\sigma_u^2 = 0$ は「グループ差が存在しない」という結論ではありません。推定値が境界にあるという状態です。$\hat\sigma_u^2$ に依存する ICC やランダム効果の予測値も、ほぼ 0 の値になります。また、分散成分に関する漸近的な近似は真のパラメータが境界の内部にあることを前提とするため、境界解ではこの前提が成り立ちません。この状態の検出と対処については GLMM タブのページを参照してください。

AIC/BIC の制限

固定効果が異なるモデルを AIC/BIC で比較できるのは、Laplace 近似で推定する組み合わせだけです。REML の射影は固定効果の構成に依存するため、REML ベースの AIC/BIC は固定効果が異なるモデル間で比較できません。MIDAS が REML で推定するのは Gaussian + identity だけなので、この制約を受けるのはこの組み合わせだけです。それ以外の組み合わせでは、Laplace 近似した周辺対数尤度を固定効果も含むすべてのパラメータについて最大化した値から AIC/BIC を導出します。この値は射影した残差ではなく応答 $y$ の周辺尤度の近似で、固定効果の構成によって尤度の定義が変わらないため、固定効果が異なるモデルどうしを比較できます。Laplace 近似の誤差は固定効果の構成によっても変わるので、AIC/BIC の差には近似の誤差の差も含まれます。また、この比較は同じ観測で当てはめたモデルどうしに限られます。MIDAS は選択した列のどれかが欠損している行を除外するので、欠損を含む説明変数を足し引きすると使われる観測が変わり、対数尤度を比べられなくなります。

ただし、どちらの経路でも比較できるのは同一の分布族とリンクの中だけです。分布族やリンクが変わると、AIC/BIC の差にモデルの当てはまり以外の要因が混ざります。推定経路が変われば対数尤度の基準（REML か最尤か）が揃わず、分布族が変われば対数尤度の測る対象（離散分布の確率質量か連続分布の密度か）や分散パラメータの扱いが変わり、リンクが変われば Laplace 近似の誤差の出方が変わるためです。

ランダム効果の予測と縮小推定

当てはめの結果には、固定効果の係数だけでなく、グループごとの $u_j$ の値も含まれます。$u_j$ の値は、そのグループのデータだけから計算されるのではありません。$u_j$ には確率分布 $N(0, \sigma_u^2)$ を仮定しているため、データとその分布という 2 つの情報源を組み合わせて値を推測します。この組み合わせが、この節で説明する縮小推定を生みます。

固定効果との扱いの違いは、用語にも現れます。固定効果 $\beta$ は未知の定数なので「推定（estimation）」と呼び、$u_j$ はモデル内で確率変数として扱われるため「予測（prediction）」と呼び分けます。$\beta$ ではデータの有限性だけが不確実性の原因なのに対し、$u_j$ では仮定した分布も情報源として働く、という違いを反映した形式的な区別で、$u_j$ が本当にランダムに生成されたという主張ではありません。

Gaussian（LMM）では、ランダム効果の予測に BLUP（Best Linear Unbiased Predictor、最良線形不偏予測量）を使います。BLUP は、$u_j$ に対する不偏性 $E[\hat u_j] = E[u_j] = 0$ が未知の固定効果 $\beta$ のどの値でも成り立つという制約のもとで、平均二乗予測誤差 $E[(\hat u_j - u_j)^2]$ を最小化する線形予測量で、Henderson (1975) の混合モデル方程式から導かれます。この予測量は縮小推定（shrinkage estimator）の形をとります:

$$\hat{u}_j = \frac{n_j}{n_j + \sigma_e^2 / \sigma_u^2} \times \bar{r}_j$$

$n_j$ はグループ $j$ のサイズ、$\bar{r}_j$ はグループ $j$ の平均残差（固定効果で説明できない部分）です。

係数 $n_j / (n_j + \sigma_e^2/\sigma_u^2)$ は0から1の値をとり、グループサイズが大きいほど1に近づきます。つまり:

• 大きなグループ: データが十分にあるので、そのグループ固有の推定値をほぼそのまま使う
• 小さなグループ: データが少ないので、全体平均（ゼロ）に向かって引き寄せる

これは「情報の借用（borrowing strength）」とも呼ばれます。この言い回しは、Tukey が 1960 年代の選挙速報予測の仕事で縮小推定の考え方に与えた呼び名に由来します（Brillinger, 2002）。データの少ないグループは他のグループの情報を借りて推定を安定させます。グループ固有の推定値をそのまま使うと分散が大きくなりますが、全体平均だけを使うとグループの特性を無視してしまいます。BLUP はこのバイアスとバリアンスのトレードオフを最適にバランスさせます。

予測値には不確かさの評価が付きます。縮小の係数を $c_j = n_j / (n_j + \sigma_e^2/\sigma_u^2)$ と書くと、$\hat u_j$ の標準誤差は、観測データで条件付けた $u_j$ の分散 $\sigma_u^2 (1 - c_j)$ の平方根です。推定した分散成分は真の値として扱います。グループのデータが少ないほど $c_j$ は 0 に近づくため、予測値は 0 へ縮み、標準誤差は分布の標準偏差 $\sigma_u$ に近づきます。データが増えるほど標準誤差は 0 へ向かいます。

Gaussian + identity 以外の組み合わせでは、MIDAS はランダム効果の予測値として条件付きモード（conditional mode）を使います。観測された応答で条件付けたときの $u_j$ の分布において、密度が最大になる値のことです。仮定した分布 $N(0, \sigma_u^2)$ を事前分布とみなせば、事後分布の最頻値にあたります。条件付きモードは、推定した $\beta$ と $\theta$ のもとで $D_p$ を $u_j$ について最小にする値です。条件付きモードも BLUP と同じく全体平均へ縮小しますが、線形予測量ではないので BLUP ではありません。

条件付きモードの標準誤差は、条件付きモードの周りの正規近似の分散 $\phi/(C_j + 1/\theta^2)$ の平方根です。$C_j$ は Laplace 近似の式の $C_j$ と同じ量です。推定した $\beta$・$\theta$・$\phi$ は真の値として扱います。

ICC（級内相関係数）

ICC（Intraclass Correlation Coefficient）は、固定効果で説明される分を除いた残差側の分散のうち、グループ間の違いが占める割合です:

$$\text{ICC} = \frac{\sigma_u^2}{\sigma_u^2 + \sigma_e^2}$$

$\sigma_u^2$ と $\sigma_e^2$ は当てはめたモデルの分散成分です。固定効果 $X\beta$ を含むモデルでは、どちらも固定効果で説明された分散を除いた残差側の量になります。

ICC は、同一グループ内の 2 観測の相関でもあります。ランダム切片モデルの節で導いたとおり、同じグループの 2 観測の相関は $\sigma_u^2 / (\sigma_u^2 + \sigma_e^2)$ で、この定義式と同じ量です。級内相関係数という名前はこの性質を指しています。

この形の ICC を条件付き ICC（conditional ICC）と呼び、MIDAS の GLMM が計算する ICC もこれです。対になる無条件 ICC は、説明変数を含めない切片のみのモデルから求めます。条件付き ICC は固定効果の指定に依存するため、説明変数を足し引きすると値が変わり、無条件 ICC とは一般に一致しません。

Non-Gaussian の場合

Non-Gaussian 分布族では $\sigma_e^2$ をデータから直接推定できません。ICC が分散分解として成立するには、$\sigma_u^2$ と $\sigma_e^2$ が同じ尺度上の同じ量のばらつきを分け合う成分である必要があります。この条件を満たす $\sigma_e^2$ は、リンク関数と分布の仮定から理論的に導出できる場合にしか得られません。

Binomial + logit/probit には閾値モデル（Goldstein et al. 2002）があります。観測されない連続的な潜在変数 $y^* = X\beta + u + e$ が $y = \mathbf{1}(y^* > 0)$ として二値応答を生成します。$e$ の分布はリンク関数の選択から演繹的に決まります。logit リンクでは $e$ がロジスティック分布に従い分散は $\pi^2/3$、probit リンクでは $e$ が $N(0, 1)$ に従い分散は $1$ です。この枠組みでは $\sigma_u^2 + \sigma_e^2$ が潜在変数の全分散を表し、ICC は分散分解として成立します。

Binomial の ICC は潜在変数 $y^*$ の尺度上の値です。観測された二値応答どうしのグループ内相関はこれとは別の量で、二値化で情報が落ちる分だけ潜在尺度の値より小さく、全体の発生確率にも依存します。

ICC が分散分解として成立するのは Gaussian + identity（REML で $\sigma_e^2$ を直接推定）と上記 2 つの Binomial の組み合わせのみです。それ以外の family+link では ICC が定義されません:

• Poisson（全リンク）: Poisson 分布では $\operatorname{Var}(Y \mid \mu) = \mu$ であり、分散は平均から完全に規定されます。独立した残差分散パラメータが存在しません。ICC の分母に $1$ を入れる規約もありますが、$\sigma_u^2 / (\sigma_u^2 + 1)$ は $\sigma_u^2$ の単調変換に過ぎず分散分解の解釈を持ちません。
• Gamma（全リンク）: profiled な分散パラメータ $\phi$ は条件付き分散のパラメータであり、潜在尺度の残差分散ではありません。
• Gaussian + log: $\sigma_u^2$ はリンク尺度上の量で、$\sigma_e^2$ は応答尺度で推定された量です。両者を同じ尺度に揃える理論残差分散が存在しません。

デザイン効果と GLM との使い分け

グループ構造を無視した GLM で標準誤差がどれだけ過小になるかは、ICC とグループサイズから見積もれます。同じグループ内の観測は互いに似ていて独立な観測ほどの情報を持たないため、係数の推定量の分散は、同じサイズの独立なデータから推定した場合より大きくなります。この倍率がデザイン効果（DEFF）で、$\text{DEFF} = 1 + (\bar n - 1) \times \text{ICC}$（$\bar n$ はグループの平均サイズ）と表せます。この式はグループサイズが均等な場合に正確で、不均等な場合は近似値です。独立を仮定した標準誤差の計算はこの増加を織り込まないため、真の標準誤差は報告値のおよそ $\sqrt{\text{DEFF}}$ 倍です。

分散がどれだけ増えるかは説明変数によって異なります。グループ平均やグループレベルで変動する説明変数の係数ではこの目安がそのまま当てはまり、グループ内で変動する説明変数の係数では増加は小さく、グループ間の変動をほとんど持たない説明変数ではほとんど増えません。

この見積もりは GLM と混合モデルの使い分けの材料になります。ICC が小さくグループサイズも大きくなければ DEFF は 1 に近く、グループ構造を無視して GLM で分析してもほぼ同じ結果が得られます。DEFF が 1 を大きく超えるなら、混合モデルが必要です。

固定効果の推測

固定効果の信頼区間の構成は、LMM（Gaussian + identity）とそれ以外の GLMM で異なります。LMM では分散成分を推定したことによる不確実性を Kenward-Roger 法で反映し、それ以外の GLMM では標準正規近似を使います。

LMM: Kenward-Roger 法

LMM の固定効果 $\hat\beta_j$ の信頼区間は $\hat\beta_j \pm t_{1-\alpha/2,\,\nu_j} \cdot \text{SE}(\hat\beta_j)$ として構成します。素朴な標準誤差は $(X'\hat V^{-1}X)^{-1}$（$\hat V = \hat\sigma_u^2 ZZ' + \hat\sigma_e^2 I$、$Z$ は各観測のグループ所属を表す指示行列）の対角成分から取れますが、この量は分散成分 $(\sigma_u^2, \sigma_e^2)$ を推定値で置き換えたことによる追加のばらつきを含まず、標準誤差を過小に見積もります。Kenward-Roger 法はこの過小分を補正した調整共分散行列から標準誤差を取り、あわせて係数ごとの自由度 $\nu_j$ を分散成分の推定精度から導きます（Kenward & Roger, 1997）。

自由度 $\nu_j$ は、その係数がどれだけの情報量で推定されているかを表します。グループ間でしか変動しない説明変数の係数は、実効的な標本サイズがグループ数で決まるため自由度が小さくなります。グループレベルの説明変数だけを均衡した設計で用いた場合、その係数の推測はグループ平均に対する通常の回帰と一致し、自由度もその残差自由度（グループ数 − 推定するグループレベル係数の数）に一致します。一方、グループ内で変動する説明変数の係数は自由度が大きく、t 分布は標準正規に近づきます。自由度が小さい係数では t 分位点が大きくなり、信頼区間はその分広がります。

グループ数が少ないデータでは、補正後の信頼区間は正規近似の区間よりはっきり広くなります。たとえば自由度 3 の t 分位点 $t_{0.975,3} \approx 3.18$ は $z_{0.975} \approx 1.96$ の約 1.6 倍です。この広さは分散成分が少数のグループからしか推定できないという情報の少なさをそのまま表したもので、正規近似の狭い区間は名目の被覆確率を下回ります。

この構成は、応答の正規性、ランダム切片の正規性、グループ間の独立、そして共分散構造 $V = \sigma_u^2 ZZ' + \sigma_e^2 I$ が正しいことに依存します。これらはいずれも緩められない仮定です。仮定が崩れている場合、たとえばランダム傾きが必要なデータをランダム切片モデルで分析している場合には、補正の有無によらず標準誤差自体が信頼できません（ランダム傾きと交差ランダム効果を参照）。

分散成分の構造が退化していて調整を計算できない場合、MIDAS は警告を表示し、未調整の標準誤差と標準正規近似で信頼区間を構成します。

Non-Gaussian GLMM: 標準正規近似

Gaussian + identity 以外の組み合わせでは、信頼区間は $\hat\beta_j \pm z_{1-\alpha/2} \cdot \text{SE}(\hat\beta_j)$ として構成します。$\text{SE}(\hat\beta_j)$ は $(X'\hat V^{-1}X)^{-1}$ の対角成分の平方根です。$\hat V = \hat\phi(\hat W^{-1} + \hat\theta^2 ZZ')$ は GLM の IRLS の調整従属変数の近似共分散で、Poisson と Binomial では $\hat\phi = 1$ です。$\hat W$ は IRLS の重み行列を、固定効果の推定値と条件付きモードから計算した平均 $\hat\mu$ で評価したものです。$\hat\phi$ は GLM の分散パラメータ（Gaussian の deviance/(n−p)、Gamma の Pearson $\chi^2$/(n−p)）とは定義が異なるので、グループ間の分散がほぼ 0 のデータでも、GLMM と GLM の標準誤差は一致しません。この構成は、$(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ が漸近的に標準正規分布に従うことを根拠とし、分散成分の推定値を真の値として扱います。

分散成分を推定値で置き換えると追加の不確実性が生じ、$(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ の分布は標準正規より裾が重くなります。グループ数が多ければ分散成分の推定が安定し、この近似は十分に正確です。グループ数が少ない場合は信頼区間が実際より狭くなる傾向があります。$\hat\phi$ も、自由度を補正しない最尤推定値なので、観測数が少ないと $\phi$ を小さく見積もり、信頼区間を狭くする方向に働きます。

ランダム傾きと交差ランダム効果

MIDAS の GLMM が対応するのはランダム切片モデルだけです。ランダム切片が表すのはグループごとのベースラインの違いまでなので、この節では対応範囲の外にある 2 つの拡張を、ランダム切片モデルで足りるかどうかを判断する材料として説明します。

説明変数の効果もグループごとに異なる場合は、ランダム傾き（random slope）モデルが必要です:

$$g(\mu_i) = x_i'\beta + u_{0j[i]} + u_{1j[i]} x_{1i}$$

$u_{0j}$ はランダム切片、$u_{1j}$ はランダム傾きで、両者は多変量正規分布に従います。

傾きに群間変動があるデータをランダム切片モデルで済ませると、その変動は誤差項に入ります。このとき群内の共分散は実際には説明変数の値に依存しますが、ランダム切片モデルはそれを定数 $\sigma_u^2$ で置き換えます。共分散構造の誤特定が残るため、当の説明変数の係数の標準誤差は信頼できなくなります。

観測が入れ子になっていない複数のグルーピング変数に属する場合は、交差ランダム効果（crossed random effects）を使います。たとえば生徒が学校と地域の両方に属するなら、グルーピング変数ごとのランダム効果を線形予測子に並べて加えます:

$$g(\mu_i) = x_i'\beta + u_{j[i]} + v_{k[i]}, \quad u_j \sim N(0, \sigma_u^2), \quad v_k \sim N(0, \sigma_v^2)$$

$j[i]$ と $k[i]$ は観測 $i$ が属する学校と地域を表す添字で、2 つのランダム効果は互いに独立と仮定します。

ランダム切片モデルはグルーピング変数を 1 つしか持たないため、交差構造のデータでは、もう一方のグルーピング変数に由来する変動が残差に残ります。残った変動はその変数のグループ内に相関を生み、階層データと独立性の問題の節で述べたのと同じ形で、無視すれば標準誤差が過小になります。

See also

• GLMM タブ - GLMM の操作方法と結果の読み方
• GLM の基礎 - GLMM の基礎となる一般化線形モデルの理論
• OLS の基礎 - 線形モデルの数理的背景
• 統計用語集 - 用語の定義

参考文献

• Gelman, A. (2005). Analysis of variance—why it is more important than ever. The Annals of Statistics, 33(1), 1-53. https://www.jstor.org/stable/3448650
• Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1-48. https://www.jstatsoft.org/v67/i01/
• Henderson, C. R. (1975). Best linear unbiased estimation and prediction under a selection model. Biometrics, 31(2), 423-447. https://www.jstor.org/stable/2529430
• Kenward, M. G., & Roger, J. H. (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics, 53(3), 983-997. https://www.jstor.org/stable/2533558
• Brillinger, D. R. (2002). John W. Tukey: His life and professional contributions. The Annals of Statistics, 30(6), 1535-1575. https://doi.org/10.1214/aos/1043351246
• Goldstein, H., Browne, W., & Rasbash, J. (2002). Partitioning variation in multilevel models. Understanding Statistics, 1(4), 223-231. https://doi.org/10.1207/S15328031US0104_02

このページの Markdown 版もあります。