GLMM の基礎

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

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

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

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

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

この誤特定には2つの側面があります。

共分散構造の誤り: モデルは観測間の共分散をゼロと仮定しますが、同じグループ内の観測には正の相関があります。この相関を無視すると、固定効果の標準誤差と予測区間がともに過小に計算されます。95% 信頼区間は実際には 95% の被覆率を持ちません。また、グループ間とグループ内のばらつきを区別できないため、データの階層構造を定量的に把握できません。

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

混合モデル（mixed model）は、グループ内の相関を変量効果として明示的にモデル化することで、共分散構造の問題を直接的に解決します。省略変数バイアスに対しては、交絡要因を固定効果として投入する必要があります。

固定効果と変量効果

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

固定効果（fixed effects） は、研究者が推定したい体系的な関係です。たとえば「勉強時間が1時間増えるとテスト成績が何点上がるか」という効果は固定効果です。

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

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

この解釈では、固定効果/変量効果の区別は「グループがランダムに抽出されたかどうか」ではなく、「グループ間の差をどうモデル化するか」の設計判断です（Gelman, 2005）。グループ差の全体的なばらつきに関心がある場合は変量効果が適切で、特定のグループの効果を直接推定したい場合は固定効果として扱います。

ランダム切片モデル

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

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

$g$ はリンク関数、 $x_i'\beta$ は固定効果の線形予測子、 $u_{j[i]}$ は観測 $i$ が属するグループ $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$ はグループ内（個人レベル）のばらつきです。

パラメータ推定

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

REML（制限付き最尤法）

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

REML は $\beta$ を射影で除去した残差の尤度を最大化することで、この偏りを補正します。

AIC/BIC の制限

REML の射影は固定効果の構成に依存するため、REML ベースの AIC/BIC は固定効果が異なるモデル間で比較できません。MIDAS の GLMM は AIC/BIC を表示しますが、この制約が適用されます。Laplace 近似の対数尤度から導出した AIC/BIC には近似誤差が伝播するため、同一ファミリー・リンク内での比較に限定してください。

Laplace 近似と PIRLS

Gaussian ファミリーと identity リンクの組み合わせでは変量効果の積分が解析的に実行できるため REML で十分ですが、それ以外の組み合わせ（Binomial, Poisson, Gamma, および Gaussian + log）では、変量効果を積分除去する操作:

L(\beta, \sigma_u^2) = \int \prod_i f(y_i \mid \beta, u) \cdot \phi(u; 0, \sigma_u^2) \, du

に閉じた形の解がありません。 $\phi$ は正規分布の密度関数です。

Laplace 近似はこの積分を、被積分関数の最大点周辺での2次近似で置き換えます。MIDAS は Bates et al. (2015) が記述した lme4 のアルゴリズムに基づいています。Laplace 近似はグループあたりの観測数が少ない場合や、Binomial でイベントが稀な場合に精度が低下することがあります。

具体的には:

外側ループ: 相対共分散パラメータ $\theta$ を Nelder-Mead 法で最適化
内側ループ: PIRLS（Penalized IRLS）で $(\beta, u)$ を同時推定

PIRLS は GLM の IRLS にランダム効果のペナルティ項を加えたものです。ペナルティ付き逸脱度（Bates et al., 2015 の表記に従い、 $\theta$ は相対共分散パラメータで、Gaussian では $\theta = \sigma_u/\sigma_e$ 、Poisson/Binomial では $\theta = \sigma_u$ 、Gamma では $\theta = \sigma_u/\sqrt\phi$ ）:

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

第1項はデータへの適合度、第2項は $u_j$ が正規分布から大きく外れることへのペナルティです。このペナルティにより、データの少ないグループの推定値は全体平均に引き寄せられます（後述の縮小推定を参照）。

推定と予測

固定効果 $\beta$ は未知の定数なので「推定（estimation）」の対象です。一方、変量効果 $u_j$ はモデル内で確率変数として扱われるため、形式的には「予測（prediction）」の対象になります。

$u_j$ が「予測」と呼ばれるのは、モデルの数学的な構造に由来します。 $u_j$ に確率分布を仮定しているため、データだけでなくその分布からの情報も組み合わせて $u_j$ の値を推測します。 $\beta$ の推定ではデータの有限性だけが不確実性の原因ですが、 $u_j$ の予測ではモデルが仮定した分布（ $N(0, \sigma_u^2)$ ）も情報源として働きます。この2つの情報源の組み合わせが、後述する縮小推定を生み出します。

「予測」という用語は、 $u_j$ が本当にランダムに生成されたことを前提としているわけではありません。グループ差を確率分布で近似するというモデリング上の判断に伴う、形式的な区別です。

変量効果の予測に BLUP（Best Linear Unbiased Predictor）がよく使われます。BLUP は、固定効果の推定可能な線形関数が不偏であるという制約のもとで、平均二乗予測誤差 $E[(\hat u_j - u_j)^2]$ を最小化する線形予測量で、Henderson の混合モデル方程式から導かれます。

Gaussian の場合、 $u_j$ の BLUP は縮小推定（shrinkage estimator）の形をとります。Non-Gaussian の場合、変量効果は条件付きモード（事後分布の最頻値）として推定されます。これは PIRLS の反復解として得られ、Gaussian の BLUP と同様に全体平均への縮小を示しますが、厳密には BLUP ではありません。

Gaussian での BLUP の式:

\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）」とも呼ばれます。データの少ないグループは他のグループの情報を借りて推定を安定させます。グループ固有の推定値をそのまま使うと分散が大きくなりますが、全体平均だけを使うとグループの特性を無視してしまいます。BLUP はこのバイアスとバリアンスのトレードオフを最適にバランスさせます。

ICC（級内相関係数）

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

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

固定効果 $X\beta$ を含むモデルでは、分子・分母とも固定効果で説明された分散を除いた残差側の量であり、これは条件付き ICC です。

ICC = 0 ならグループ間に差がなく、すべてのばらつきは個人レベルです。ICC = 1 ならグループ内の個人差がなく、すべてのばらつきはグループ間の差です。

ICC が小さければ、グループ構造を無視して GLM で分析してもほぼ同じ結果が得られます。ただし影響の大きさは ICC だけでなくグループサイズにも依存します。設計効果 $\text{DEFF} = 1 + (\bar n - 1) \times \text{ICC}$ （ $\bar n$ はグループの平均サイズ）が分散の膨張倍率の目安になります。この式はグループサイズが均等な場合に正確で、不均等な場合は近似値です。グループ平均やグループレベルで変動する予測変数の係数では、標準誤差の目安が $\sqrt{\text{DEFF}}$ 倍に膨らみます。グループ内で変動する予測変数の係数では膨張の度合いが異なります。ICC が大きい場合は混合モデルを使う必要があります。

Non-Gaussian の場合

Non-Gaussian ファミリーでは $\sigma_e^2$ をデータから直接推定できません。ICC を分散分解として成立させるには、 $\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 + Logit	$\pi^2/3 \approx 3.29$	ロジスティック分布（閾値モデル）
Binomial + Probit	$1$	標準正規分布（閾値モデル）

Binomial の ICC は潜在尺度上の値であり、観測された二値応答の相関とは一致しません。サンプルサイズ設計に使う場合はこの違いに注意してください。

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

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

MIDAS の実装

MIDAS の GLMM は分布ファミリーとリンク関数の組み合わせに応じて異なるアルゴリズムを使います。

Gaussian + identity（LMM）: Profile REML。分散成分の比 $\theta^2 = \sigma_u^2/\sigma_e^2$ を Nelder-Mead 法で数値最適化し、各 $\theta$ に対して $\beta$ と $\sigma_e^2$ を閉じた形で求めます。

それ以外の組み合わせ（Binomial, Poisson, Gamma, および Gaussian + log）: Bates et al. (2015) に基づく Laplace 近似。外側ループで相対共分散パラメータ $\theta$ を Nelder-Mead 法で最適化し、内側ループで PIRLS により $(\beta, u)$ を同時推定します。分散パラメータ $\phi$ は Poisson/Binomial では1に固定します。Gamma と Gaussian + log では、最適化中の尤度評価に条件付き逸脱度から求めた $\phi = \sum_i d(y_i, \mu_i)/n$ を使い、結果として表示する $\phi$ は Pearson $\chi^2/(n-p)$ で推定します。

Laplace 近似の経路では、固定効果の初期値に GLM（ランダム効果なし）の推定結果を使います。説明変数は内部的にスケーリング（標準化）され、推定後に元のスケールに逆変換されます。

固定効果の推測と正規近似

固定効果 $\hat\beta_j$ の信頼区間は Wald 法に基づき、 $\hat\beta_j \pm z_{1-\alpha/2} \cdot \text{SE}(\hat\beta_j)$ として構成します。この構成は、 $(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ が近似的に標準正規分布に従うことを根拠とします。

通常の回帰（OLS）では誤差分散 $\sigma^2$ を推定するため $(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ は t 分布に従い、自由度は $n - p$ です。混合モデルでは分散成分が複数あり ( $\sigma_u^2$ , $\sigma_e^2$ )、対応する自由度が一意に定まりません。そのため分散成分の推定値を真の値として扱い、標準正規分布で近似します。

LMM では $u_j$ と $\varepsilon_i$ に正規分布を仮定しているため、分散成分が既知なら $(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ は正規分布に正確に従います。Non-Gaussian GLMM では正規分布は漸近的な近似です。いずれの場合も、分散成分を推定値で置き換えると追加の不確実性が生じ、 $(\hat\beta_j - \beta_j) / \text{SE}(\hat\beta_j)$ の分布は標準正規より裾が重くなります。

グループ数が多ければ分散成分の推定が安定し、この近似は十分に正確です。グループ数が少ない場合は信頼区間が実際より狭くなる傾向があります。

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

ランダム切片モデルはグループごとのベースラインの違いだけをモデル化します。説明変数の効果もグループごとに異なる場合は、ランダム傾き（random slope）モデルが必要です:

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

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

また、観測が複数のグルーピング変数に属する場合（たとえば生徒が学校と地域の両方に属する）、交差ランダム効果（crossed random effects）が使えます。

MIDAS は現在ランダム切片モデルのみ対応しています。

参考文献

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/
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 版もあります。