階層回帰と混合効果モデル

🎓 レベル：発展　|　重要度：B（標準）

📎 前提：ベイズGLM　|　関連：階層モデルの構造・ベイズ線形回帰

要点（BLUF）

階層回帰（混合効果モデル）は、群ごとに切片・傾きが変わる回帰を、それらの係数が共通の親分布から生まれたとみなして結びます。第5章の部分プーリングを回帰に持ち込んだもの。
変動切片モデル $y_{ij}=\alpha_j+\beta x_{ij}+\varepsilon,\ \alpha_j\sim\mathcal N(\mu,\tau^2)$ 。群切片が部分プーリングで全体平均へ収縮し、データの少ない群ほど強く収縮します。
ギブスサンプリング（第4章）で全パラメータを同時推定。共通の $\beta$ （固定効果）＋群ごとの $\alpha_j$ （変動効果）＝混合効果モデルです。

1. 群構造のある回帰

店舗ごと・学校ごと・被験者ごとに繰り返し測定がある回帰では、群によって水準（切片）や効果（傾き）が違います。素朴には群ごとに別々の回帰（非プーリング）か、群を無視した1本の回帰（完全プーリング）ですが、階層モデルの構造と同じく部分プーリングが両者の良いとこ取りをします。

最もよく使う変動切片モデルは、傾き $\beta$ は全群共通、切片 $\alpha_j$ だけ群ごとに変わるとし、その $\alpha_j$ を親分布で結びます。

y_{ij}\sim\mathcal N(\alpha_j+\beta x_{ij},\ \sigma^2),\qquad \alpha_j\sim\mathcal N(\mu,\tau^2)

$\beta$ のように全群共通の係数を固定効果、 $\alpha_j$ のように群ごとに変わる係数を変動効果（ランダム効果）と呼び、両方を持つので混合効果モデルです。

2. ギブスで全パラメータを同時推定する

このモデルの完全条件付きは、すべて第2章の共役（正規–正規・逆ガンマ）になり、ギブス（ギブスサンプリング）が回せます。 $\alpha_j\mid\text{rest}$ は群 $j$ の残差 $y_{ij}-\beta x_{ij}$ に対する正規–正規で、精度の加重平均で全体平均 $\mu$ へ収縮します。

import numpy as np
from scipy import stats

# 変動切片モデル：y_ij = α_j + β x_ij + ε、α_j~N(μ,τ²)
rng = np.random.default_rng(3)
J = 8
n_j = np.array([4, 5, 6, 8, 12, 20, 30, 40])           # 群サイズ（小～大）
beta_true, mu_true, tau_true, sigma = 1.5, 2.0, 1.5, 1.0
alpha_true = rng.normal(mu_true, tau_true, J)
groups = []
for j in range(J):
    x = rng.uniform(-2, 2, n_j[j])
    y = alpha_true[j] + beta_true*x + rng.normal(0, sigma, n_j[j])
    groups.append((x, y))

# 非プーリング：各群で切片を単純推定（共通βを仮定し群平均残差）
np_alpha = np.array([np.mean(y - 1.5*x) for (x, y) in groups])

# ギブス：α_j, β, μ, τ², σ² を同時サンプリング
def gibbs(N=40000, burn=8000, seed=0):
    rng = np.random.default_rng(seed)
    alpha = np_alpha.copy(); beta = 1.0; mu = alpha.mean(); tau2 = s2 = 1.0
    allx = np.concatenate([g[0] for g in groups]); ally = np.concatenate([g[1] for g in groups])
    gidx = np.concatenate([[j]*n_j[j] for j in range(J)])
    A = np.zeros((N, J))
    for it in range(N):
        r = ally - alpha[gidx]                                  # β | rest
        vb = 1/(allx@allx/s2 + 1e-6); beta = rng.normal(vb*(allx@r/s2), np.sqrt(vb))
        for j in range(J):                                      # α_j | rest（収縮）
            x, y = groups[j]; rj = (y - beta*x)
            prec = n_j[j]/s2 + 1/tau2
            alpha[j] = rng.normal((rj.sum()/s2 + mu/tau2)/prec, np.sqrt(1/prec))
        mu = rng.normal(alpha.mean(), np.sqrt(tau2/J))          # μ | rest
        tau2 = stats.invgamma(a=1+J/2, scale=1+0.5*np.sum((alpha-mu)**2)).rvs(random_state=rng)
        resid = ally - alpha[gidx] - beta*allx
        s2 = stats.invgamma(a=1+len(ally)/2, scale=1+0.5*resid@resid).rvs(random_state=rng)
        A[it] = alpha
    return A[burn:]

pp_alpha = gibbs().mean(0)
print(f"{'群':<3}{'n_j':>5}{'真α_j':>8}{'非プーリング':>12}{'部分プーリング':>14}{'収縮量':>9}")
for j in range(J):
    print(f"{j+1:<3}{n_j[j]:>5}{alpha_true[j]:>8.2f}{np_alpha[j]:>12.2f}{pp_alpha[j]:>14.2f}{np_alpha[j]-pp_alpha[j]:>+9.2f}")

出力：

群    n_j    真α_j      非プーリング       部分プーリング      収縮量
1      4    5.06        4.46          4.20    +0.26
2      5   -1.83       -1.44         -1.21    -0.23
3      6    2.63        2.48          2.43    +0.05
4      8    1.15        1.51          1.46    +0.04
5     12    1.32        1.25          1.25    +0.00
6     20    1.68        1.85          1.85    +0.01
7     30   -1.03       -0.90         -0.86    -0.04
8     40    1.65        1.56          1.55    +0.02

出力の意味：群切片 $\alpha_j$ が部分プーリングで全体へ収縮します。データの少ない群1（ $n=4$ ）・群2（ $n=5$ ）は収縮量が $\pm0.25$ 前後と大きく、全体平均の方へ引き戻されます。一方、データの多い群（ $n\ge20$ ）は収縮量がほぼ $0$ で、自前の推定をほぼ保ちます。階層モデルの構造の収縮が、回帰の切片に対してそのまま働いているわけです。全群共通の傾き $\beta$ （固定効果）と群ごとの切片 $\alpha_j$ （変動効果）を、ギブスが一度に推定しています。

3. 固定効果・変動効果・部分プーリング

固定効果（共通 $\beta$ ）：全群に共通の効果。プールして強く推定。
変動効果（群ごと $\alpha_j$ ）：群差。親分布 $\mathcal N(\mu,\tau^2)$ で結び、部分プーリングで収縮。
$\tau$ （群間ばらつき）の大小が収縮の強さを決める（ $\tau\to0$ で完全プーリング、 $\tau\to\infty$ で非プーリング）——階層モデルの構造と同じ構図。

切片だけでなく傾きも群ごとに変える（変動傾きモデル $y_{ij}=\alpha_j+\beta_j x_{ij}$ 、 $(\alpha_j,\beta_j)$ を多変量正規で結ぶ）と表現力が上がります。ただし階層回帰の事後は漏斗型になりやすく、中心化のままだとサンプラーが詰まるので、非中心化が要ります（階層モデルの実例と再パラメータ化）。非ガウスの応答（二値・カウント）なら階層 GLM＝ベイズGLM にランダム効果を足した形で、実務では PyMC などで NUTS が定番（要最新確認）。

⚠️ よくある誤解

「群ごとに別回帰（非プーリング）が一番正確」ではない：小群は推定が暴れます。部分プーリングが安定（収縮の数理）。
「混合効果は固定効果だけ／変動効果だけ」ではない：共通の固定効果と群ごとの変動効果を両方持つのが混合効果モデル。
「収縮は全群一律」ではない：データの少ない群ほど強く収縮（§2）。 $\tau$ と群サイズで決まります。
「階層回帰はギブスで素直に回る」とは限らない：漏斗型で詰まることがあり、非中心化が要る（階層モデルの実例と再パラメータ化）。

まとめ（Phase 7）

第7章では、回帰をベイズで組み立てました——係数の事後と予測帯が出るベイズ線形回帰（ベイズ線形回帰）、正則化＝MAP の対応（正則化との関係）、リンク関数で共役が崩れる GLM をラプラス近似で解く（ベイズGLM）、そして群構造を収縮で扱う階層回帰（本ノート）。第2章の共役・第4章の MCMC・第5章の階層が、回帰という応用で一つに合流しました。次章では、組んだモデルが良いかどうかを評価・比較する道具へ進みます。