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 は固定効果が異なるモデル間で比較できません。固定効果の構成が異なるモデルを比較するには、ML で再推定した AIC/BIC を使う必要があります。

Laplace 近似と PIRLS

Gaussian ファミリーでは変量効果の積分が解析的に実行できるため REML で十分ですが、Non-Gaussian ファミリー（Binomial, Poisson, Gamma など）では、変量効果を積分除去する操作:

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}

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$ が直接定義されないため、潜在変数の分散を代用します。Binomial + Logit リンクでは、観測される二値応答の背後に連続的な潜在変数があると考えます。Logit リンクの場合、この潜在変数はロジスティック分布に従い、その分散 $\pi^2/3$ を残差分散の代用とします。

ファミリー + リンク	潜在尺度の残差分散
Binomial + Logit	$\pi^2/3 \approx 3.29$
Poisson + Log	$1$ （対数正規-Poisson 混合モデルにおける近似。 $\mu$ が小さい場合、精度は低下します）
Gamma + Log	$\phi$ （推定された分散パラメータ）

MIDAS の実装

MIDAS の GLMM は分布ファミリーに応じて異なるアルゴリズムを使います。

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

Non-Gaussian（Binomial, Poisson, Gamma 等）: Bates et al. (2015) に基づく Laplace 近似。外側ループで相対共分散パラメータ $\theta$ を Nelder-Mead 法で最適化し、内側ループで PIRLS により $(\beta, u)$ を同時推定します。分散パラメータ $\phi$ は Poisson/Binomial では1に固定、Gamma では profiled deviance（ $D_p/n$ ）から推定します。

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

固定効果の検定と正規近似

固定効果 $\hat\beta_j$ の検定には Wald 検定を使います。検定統計量は $z = \hat\beta_j / \text{SE}(\hat\beta_j)$ で、p 値と信頼区間は標準正規分布から計算します。

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

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