全国の学校から生徒を集めて学力テストを実施したとします。生徒の成績は個人の能力だけでなく、所属する学校にも依存するでしょう。「学校」という要因は、調査対象の特定の学校だけでなく、全国のすべての学校を代表する母集団からのランダムなサンプルと見なすのが自然です。このような「ランダムにサンプリングされたグループ」の効果をモデルに組み込むにはどうすればよいでしょうか。
あるいは、臨床試験で同じ患者から複数回の測定を行う場合を考えてみてください。患者ごとにベースラインの値が異なり、投薬効果も個人差があるかもしれません。各測定は同じ患者に「入れ子」になっている（nested）ため、データの独立性の仮定が崩れます。通常の回帰分析を適用すると、標準誤差が過小推定され、偽陽性の原因になります。
混合効果モデル（Mixed Effects Model）は、固定効果（fixed effects）とランダム効果（random effects）の両方を含む統計モデルです。固定効果は研究者が関心のある体系的な要因（性別、処理群など）を表し、ランダム効果はランダムにサンプリングされたグループ（学校、被験者、地域など）に由来する変動を表します。
混合効果モデルを理解すると、以下のような応用が開けます。
- 教育学: 学校（クラスター）内の生徒のネスト構造を考慮した分析
- 臨床試験: 反復測定データ（同じ被験者からの複数回の測定）の分析
- 生態学: 複数のサイトから集めたデータの分析。サイト間の変動をランダム効果として扱う
- 言語学: 被験者と刺激の両方をランダム効果として扱う交差ランダム効果モデル
- スポーツ科学: 選手のパフォーマンスの個人差を考慮した分析
本記事では、混合効果モデルの数学的定式化を行い、固定効果とランダム効果の違いを理論的に明確にします。パラメータの推定法（MLとREML）を解説し、ランダム効果の予測（BLUP）と縮小推定の仕組みを導出します。Pythonで実装して確認します。
本記事の内容
- 固定効果とランダム効果の違いと使い分け
- 線形混合効果モデルの行列表現
- 周辺分布と条件付き分布
- 級内相関係数（ICC）の意味
- 最尤推定（ML）と制限付き最尤推定（REML）
- ランダム効果の予測（BLUP）と縮小推定
- Pythonでの実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
固定効果とランダム効果
直感的な理解
固定効果: 研究者が関心のある、特定の水準を持つ要因です。たとえば「薬の種類（薬A, 薬B, プラセボ）」は固定効果です。これらの水準は事前に決まっており、推論の対象そのものです。「薬Aの効果は何か？」「薬Aと薬Bの差はいくらか？」という問いに答えることが研究の目的です。
ランダム効果: 母集団からランダムにサンプリングされた、研究者が関心の対象外の要因です。たとえば「被験者」は通常ランダム効果です。特定の被験者IDのパラメータそのものには関心がなく、被験者間の変動（分散成分）に関心があります。「被験者間のばらつきはどの程度あるか？」が問いの焦点です。
この区別を料理のアナロジーで考えてみましょう。3種類の調理法（固定効果）の効果を調べるために、10人のシェフ（ランダム効果）に料理を作ってもらう実験を行ったとします。3種類の調理法は研究の対象であり、「調理法Aは調理法Bより美味しい」という結論を導きたいのです。一方、10人のシェフは母集団のランダムサンプルであり、個々のシェフの腕前そのものには関心がなく、「シェフ間のばらつきはどの程度か」を知りたいのです。
使い分けの基準
| 基準 | 固定効果 | ランダム効果 |
|---|---|---|
| 水準の性質 | 特定の水準に関心がある | 母集団からのランダムサンプル |
| 推論の対象 | 各水準の効果の値 | 水準間の分散 |
| 水準数 | 通常少ない（2〜数個） | 通常多い（10以上） |
| 反復可能性 | 同じ水準で再実験可能 | 異なるサンプルが得られる |
| 一般化の範囲 | 含まれる水準のみ | 母集団全体に一般化 |
なぜ固定効果だけでは不十分か
すべてのグループ変数を固定効果として扱うと、以下の問題が生じます。
- パラメータ数の爆発: 100の学校をダミー変数で表すと99個の追加パラメータが必要になり、自由度が大幅に減少します
- 小グループの不安定な推定: データが少ないグループの効果が不安定に推定されます。3人しかいないグループの平均は信頼できません
- 新しいグループへの予測不能: 訓練データに含まれないグループ（新しい学校）の予測ができません
- 独立性の仮定の違反: 同じグループ内のデータは共通の要因を共有するため、独立ではありません。これを無視すると標準誤差が過小推定されます
ランダム効果として扱うと、グループの効果は正規分布 $N(0, \sigma_u^2)$ からのサンプルとしてモデル化され、推定すべきパラメータは分散 $\sigma_u^2$ だけです。データが少ないグループの効果は全体平均に向かって縮小推定（shrinkage）され、安定した推定が得られます。この「部分的プーリング（partial pooling）」の概念は、混合効果モデルの最も重要な特徴の一つです。
モデルの数学的定式化に進みましょう。
線形混合効果モデルの定式化
行列表現
線形混合効果モデルは次のように書けます。
$$ \begin{equation} \bm{y} = \bm{X}\bm{\beta} + \bm{Z}\bm{u} + \bm{\epsilon} \end{equation} $$
各項の意味は以下の通りです。
- $\bm{y} \in \mathbb{R}^n$: 応答変数ベクトル（全データの目的変数）
- $\bm{X} \in \mathbb{R}^{n \times p}$: 固定効果の計画行列（切片、処理群ダミー、連続変数など）
- $\bm{\beta} \in \mathbb{R}^p$: 固定効果のパラメータベクトル（未知の定数）
- $\bm{Z} \in \mathbb{R}^{n \times q}$: ランダム効果の計画行列（グループの所属を示すインジケータ行列）
- $\bm{u} \in \mathbb{R}^q$: ランダム効果ベクトル（未知の確率変数）
- $\bm{\epsilon} \in \mathbb{R}^n$: 誤差ベクトル
$\bm{\beta}$ と $\bm{u}$ の根本的な違いに注意してください。$\bm{\beta}$ は未知の定数であり、「推定」の対象です。一方、$\bm{u}$ は確率変数であり、実現値の「予測」の対象です。
分布の仮定
ランダム効果と誤差項の分布を仮定します。
$$ \begin{equation} \bm{u} \sim N(\bm{0}, \bm{G}), \quad \bm{\epsilon} \sim N(\bm{0}, \bm{R}) \end{equation} $$
$$ \begin{equation} \text{Cov}(\bm{u}, \bm{\epsilon}) = \bm{O} \end{equation} $$
ここで $\bm{G}$ はランダム効果の共分散行列、$\bm{R}$ は誤差の共分散行列です。
最も単純なランダム切片モデルでは次のようになります。
- $\bm{G} = \sigma_u^2\bm{I}_q$（各グループのランダム効果は独立、同じ分散）
- $\bm{R} = \sigma^2\bm{I}_n$（各観測の誤差は独立、同じ分散）
ランダム傾きも含める場合は、$\bm{G}$ がブロック対角行列になります。
周辺分布の導出
ランダム効果 $\bm{u}$ を積分消去すると、$\bm{y}$ の周辺分布が得られます。
$$ \bm{y} = \bm{X}\bm{\beta} + \bm{Z}\bm{u} + \bm{\epsilon} $$
$\bm{Z}\bm{u} + \bm{\epsilon}$ は正規分布の和なので正規分布に従います。その期待値と共分散を計算します。
期待値は次のようになります。
$$ \begin{equation} E[\bm{y}] = \bm{X}\bm{\beta} + \bm{Z}E[\bm{u}] + E[\bm{\epsilon}] = \bm{X}\bm{\beta} \end{equation} $$
共分散は次のようになります。$\bm{u}$ と $\bm{\epsilon}$ が独立なので、
$$ \begin{equation} \text{Cov}(\bm{y}) = \bm{V} = \bm{Z}\bm{G}\bm{Z}^\top + \bm{R} \end{equation} $$
したがって $\bm{y} \sim N(\bm{X}\bm{\beta}, \bm{V})$ です。
級内相関係数（ICC）
共分散行列 $\bm{V} = \bm{Z}\bm{G}\bm{Z}^\top + \bm{R}$ の構造を、ランダム切片モデルで具体的に見てみましょう。
同じグループ $j$ に属する2つの観測 $y_{ij}$ と $y_{i’j}$（$i \neq i’$）の共分散は次のようになります。
$$ \text{Cov}(y_{ij}, y_{i’j}) = \text{Cov}(u_j + \epsilon_{ij}, u_j + \epsilon_{i’j}) = \text{Var}(u_j) = \sigma_u^2 $$
同じグループ内のデータは共通のランダム効果 $u_j$ を共有するため、正の共分散（相関）を持ちます。
級内相関係数（ICC: Intraclass Correlation Coefficient）はこの相関の強さを定量化します。
$$ \begin{equation} \text{ICC} = \rho = \frac{\sigma_u^2}{\sigma_u^2 + \sigma^2} \end{equation} $$
ICCは0から1の値を取り、以下のように解釈されます。
- $\text{ICC} = 0$: グループ間の変動がゼロ。グループ構造を無視してよい
- $\text{ICC} = 1$: 個人内変動がゼロ。同じグループのデータは完全に同一
- $\text{ICC} = 0.1$: 全変動の10%がグループ間、90%が個人内
一般に $\text{ICC} > 0.05$ であれば混合効果モデルの使用が推奨されます。通常の回帰分析ではグループ内の相関を無視するため、固定効果の標準誤差が過小推定され、有意性の判定が不正確になります。
パラメータの推定方法に進みましょう。
パラメータの推定
固定効果の推定（GLS）
周辺分布 $\bm{y} \sim N(\bm{X}\bm{\beta}, \bm{V})$ のもとで、$\bm{\beta}$ の一般化最小二乗推定量は次のようになります。
$$ \begin{equation} \hat{\bm{\beta}} = (\bm{X}^\top\bm{V}^{-1}\bm{X})^{-1}\bm{X}^\top\bm{V}^{-1}\bm{y} \end{equation} $$
$\bm{V}$ が既知であればこれで計算できますが、$\bm{V}$ は分散成分 $\sigma_u^2$ と $\sigma^2$ を含むため、これらも推定する必要があります。分散成分と固定効果を交互に推定する反復手続きが必要です。
最尤推定（ML）
分散成分の最尤推定は、周辺対数尤度を最大化することで行います。
$$ \begin{equation} \ell(\bm{\beta}, \bm{V}) = -\frac{n}{2}\ln(2\pi) – \frac{1}{2}\ln|\bm{V}| – \frac{1}{2}(\bm{y} – \bm{X}\bm{\beta})^\top\bm{V}^{-1}(\bm{y} – \bm{X}\bm{\beta}) \end{equation} $$
$\bm{\beta}$ を $\bm{V}$ の関数として置き換えたプロファイル尤度を分散成分について最適化します。具体的には、$\bm{V}$ が与えられたときの $\bm{\beta}$ の最適解 $\hat{\bm{\beta}}(\bm{V})$ を代入し、分散成分だけの関数にした上で最適化を行います。
制限付き最尤推定（REML）
MLは分散成分を過小推定する傾向があります。これは正規分布の標本分散 $s^2 = \sum(x_i – \bar{x})^2/n$ が母分散の不偏推定量でない（$n-1$ で割るべき）のと同じ原理です。固定効果の推定による自由度の損失を考慮しないためです。
REML（Restricted Maximum Likelihood）は、$\bm{\beta}$ に関する情報を消去した制限付き尤度を最大化することで、不偏に近い分散成分の推定量を得ます。
$$ \begin{equation} \ell_R(\bm{V}) = -\frac{n-p}{2}\ln(2\pi) – \frac{1}{2}\ln|\bm{V}| – \frac{1}{2}\ln|\bm{X}^\top\bm{V}^{-1}\bm{X}| – \frac{1}{2}\bm{r}^\top\bm{V}^{-1}\bm{r} \end{equation} $$
ここで $\bm{r} = \bm{y} – \bm{X}\hat{\bm{\beta}}$ です。
MLとREMLの違いは追加項 $-\frac{1}{2}\ln|\bm{X}^\top\bm{V}^{-1}\bm{X}|$ にあります。この項が固定効果の推定による自由度の損失を補正しています。
実用上の使い分けとして、分散成分の推定にはREMLが推奨されます。一方、固定効果の異なるモデルを尤度比検定で比較する場合にはMLを使用します（REMLは固定効果が異なるモデル間の比較に適さないため）。
ランダム効果の予測（BLUP）
ランダム効果 $\bm{u}$ は確率変数なので「推定」ではなく「予測」すると言います。最良線形不偏予測量（BLUP: Best Linear Unbiased Prediction）は、データ $\bm{y}$ が与えられたもとでの $\bm{u}$ の条件付き期待値として得られます。
多変量正規分布の条件付き分布の公式を適用します。$\bm{u}$ と $\bm{y}$ の同時分布を考えると次のようになります。
$$ \begin{pmatrix}\bm{u}\\\bm{y}\end{pmatrix} \sim N\left(\begin{pmatrix}\bm{0}\\\bm{X}\bm{\beta}\end{pmatrix}, \begin{pmatrix}\bm{G} & \bm{G}\bm{Z}^\top\\\bm{Z}\bm{G} & \bm{V}\end{pmatrix}\right) $$
条件付き期待値の公式から、BLUPは次のようになります。
$$ \begin{equation} \hat{\bm{u}} = E[\bm{u} | \bm{y}] = \bm{G}\bm{Z}^\top\bm{V}^{-1}(\bm{y} – \bm{X}\hat{\bm{\beta}}) \end{equation} $$
縮小推定の仕組み
BLUPの最も重要な性質は縮小推定です。ランダム切片モデルでグループ $j$ のBLUPを具体的に計算すると次のようになります。
$$ \begin{equation} \hat{u}_j = \frac{\sigma_u^2}{\sigma_u^2 + \sigma^2/n_j} \cdot (\bar{y}_j – \bm{\bar{x}}_j^\top\hat{\bm{\beta}}) \end{equation} $$
因子 $\gamma_j = \frac{\sigma_u^2}{\sigma_u^2 + \sigma^2/n_j}$ は縮小因子と呼ばれ、0から1の間の値を取ります。
- $n_j$ が大きい（データが多い）グループ: $\gamma_j \to 1$。グループ固有のデータを信頼し、ほぼ縮小しない
- $n_j$ が小さい（データが少ない）グループ: $\gamma_j \to 0$。グループ固有のデータを信頼せず、全体平均に強く縮小
この「データの量に応じて情報源を重み付けする」性質は、ベイズ推定の事後平均と完全に一致します。混合効果モデルのBLUPは、グループ効果に正規事前分布 $N(0, \sigma_u^2)$ を仮定したベイズ推定の経験ベイズ版と解釈できます。
Pythonで実装して確認しましょう。
Pythonによる実装と可視化
ランダム切片モデルの実装
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar
np.random.seed(42)
# --- データ生成 ---
n_groups = 20 # グループ（学校）の数
n_per_group = np.random.randint(5, 30, n_groups) # グループサイズ（不均一）
n_total = n_per_group.sum()
sigma_u_true = 2.0 # ランダム効果の標準偏差
sigma_true = 1.5 # 誤差の標準偏差
beta_true = np.array([5.0, 1.5]) # 固定効果 [切片, 傾き]
# ランダム効果（グループごとの切片のずれ）
u_true = np.random.randn(n_groups) * sigma_u_true
# データ生成
x_list, y_list, group_list = [], [], []
for j in range(n_groups):
x_j = np.random.uniform(0, 5, n_per_group[j])
y_j = beta_true[0] + u_true[j] + beta_true[1] * x_j + sigma_true * np.random.randn(n_per_group[j])
x_list.append(x_j)
y_list.append(y_j)
group_list.append(np.full(n_per_group[j], j))
x = np.concatenate(x_list)
y = np.concatenate(y_list)
group = np.concatenate(group_list).astype(int)
# --- 行列の構成 ---
X = np.column_stack([np.ones(n_total), x])
Z = np.zeros((n_total, n_groups))
for i in range(n_total):
Z[i, group[i]] = 1
# --- REML推定（1次元最適化） ---
def neg_reml(log_ratio):
ratio = np.exp(log_ratio) # sigma_u^2 / sigma^2
V = Z @ (ratio * np.eye(n_groups)) @ Z.T + np.eye(n_total)
try:
L = np.linalg.cholesky(V)
except np.linalg.LinAlgError:
return 1e10
logdet_V = 2 * np.sum(np.log(np.diag(L)))
V_inv_y = np.linalg.solve(V, y)
V_inv_X = np.linalg.solve(V, X)
XtVX = X.T @ V_inv_X
beta_hat = np.linalg.solve(XtVX, X.T @ V_inv_y)
resid = y - X @ beta_hat
V_inv_resid = np.linalg.solve(V, resid)
sigma2_hat = resid @ V_inv_resid / (n_total - X.shape[1])
reml = (n_total - X.shape[1]) * np.log(sigma2_hat) + logdet_V + np.log(np.linalg.det(XtVX))
return reml
result = minimize_scalar(neg_reml, bounds=(-5, 5), method="bounded")
ratio_hat = np.exp(result.x)
# 最終推定
V_hat = Z @ (ratio_hat * np.eye(n_groups)) @ Z.T + np.eye(n_total)
V_inv = np.linalg.inv(V_hat)
beta_hat = np.linalg.solve(X.T @ V_inv @ X, X.T @ V_inv @ y)
resid = y - X @ beta_hat
sigma2_hat = resid @ V_inv @ resid / (n_total - X.shape[1])
sigma_u2_hat = ratio_hat * sigma2_hat
# 簡易的なBLUP計算
u_hat_simple = np.zeros(n_groups)
for j in range(n_groups):
mask = group == j
n_j = mask.sum()
y_bar_j = y[mask].mean()
x_bar_j = x[mask].mean()
resid_mean = y_bar_j - beta_hat[0] - beta_hat[1] * x_bar_j
shrinkage = sigma_u2_hat / (sigma_u2_hat + sigma2_hat / n_j)
u_hat_simple[j] = shrinkage * resid_mean
icc = sigma_u2_hat / (sigma_u2_hat + sigma2_hat)
print("=== 混合効果モデルの結果 ===")
print(f"固定効果: beta = {beta_hat} (真値: {beta_true})")
print(f"sigma^2 (誤差) = {sigma2_hat:.4f} (真値: {sigma_true**2:.4f})")
print(f"sigma_u^2 (ランダム効果) = {sigma_u2_hat:.4f} (真値: {sigma_u_true**2:.4f})")
print(f"ICC = {icc:.4f} (真値: {sigma_u_true**2/(sigma_u_true**2+sigma_true**2):.4f})")
分散成分の推定値が真の値に近いことが確認できます。ICCは約0.64であり、全変動の64%がグループ間変動であることを示しています。これは混合効果モデルの使用が強く推奨される状況です。
可視化
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar
np.random.seed(42)
n_groups = 20
n_per_group = np.random.randint(5, 30, n_groups)
n_total = n_per_group.sum()
sigma_u_true = 2.0; sigma_true = 1.5
beta_true = np.array([5.0, 1.5])
u_true = np.random.randn(n_groups) * sigma_u_true
x_list, y_list, group_list = [], [], []
for j in range(n_groups):
x_j = np.random.uniform(0, 5, n_per_group[j])
y_j = beta_true[0] + u_true[j] + beta_true[1]*x_j + sigma_true*np.random.randn(n_per_group[j])
x_list.append(x_j); y_list.append(y_j); group_list.append(np.full(n_per_group[j], j))
x = np.concatenate(x_list); y = np.concatenate(y_list); group = np.concatenate(group_list).astype(int)
X = np.column_stack([np.ones(n_total), x])
Z = np.zeros((n_total, n_groups))
for i in range(n_total): Z[i, group[i]] = 1
def neg_reml(log_ratio):
ratio = np.exp(log_ratio)
V = Z @ (ratio*np.eye(n_groups)) @ Z.T + np.eye(n_total)
try: L = np.linalg.cholesky(V)
except: return 1e10
logdet_V = 2*np.sum(np.log(np.diag(L)))
Vi_y = np.linalg.solve(V, y); Vi_X = np.linalg.solve(V, X)
XtVX = X.T @ Vi_X; bh = np.linalg.solve(XtVX, X.T @ Vi_y)
r = y - X @ bh; Vi_r = np.linalg.solve(V, r)
s2 = r @ Vi_r / (n_total - 2)
return (n_total-2)*np.log(s2) + logdet_V + np.log(np.linalg.det(XtVX))
res = minimize_scalar(neg_reml, bounds=(-5,5), method="bounded")
ratio_hat = np.exp(res.x)
V_hat = Z @ (ratio_hat*np.eye(n_groups)) @ Z.T + np.eye(n_total)
Vi = np.linalg.inv(V_hat)
beta_hat = np.linalg.solve(X.T @ Vi @ X, X.T @ Vi @ y)
r = y - X @ beta_hat; sigma2_hat = r @ Vi @ r / (n_total-2)
sigma_u2_hat = ratio_hat * sigma2_hat
u_hat = np.zeros(n_groups)
u_ols = np.zeros(n_groups)
for j in range(n_groups):
mask = group == j; n_j = mask.sum()
resid_mean = y[mask].mean() - beta_hat[0] - beta_hat[1]*x[mask].mean()
shrinkage = sigma_u2_hat / (sigma_u2_hat + sigma2_hat/n_j)
u_hat[j] = shrinkage * resid_mean
Xj = np.column_stack([np.ones(n_j), x[mask]])
bj = np.linalg.lstsq(Xj, y[mask], rcond=None)[0]
u_ols[j] = bj[0] - beta_hat[0]
fig, axes = plt.subplots(1, 3, figsize=(17, 5.5))
colors = plt.cm.tab20(np.linspace(0, 1, n_groups))
# (a) グループごとの回帰直線
ax = axes[0]
for j in range(min(n_groups, 10)):
mask = group == j
ax.scatter(x[mask], y[mask], alpha=0.3, s=10, color=colors[j])
x_line = np.linspace(0, 5, 50)
y_line = (beta_hat[0] + u_hat[j]) + beta_hat[1] * x_line
ax.plot(x_line, y_line, "-", color=colors[j], linewidth=0.8, alpha=0.7)
x_line = np.linspace(0, 5, 50)
ax.plot(x_line, beta_hat[0] + beta_hat[1]*x_line, "k-", linewidth=2.5, label="Fixed effect")
ax.set_xlabel("$x$", fontsize=11); ax.set_ylabel("$y$", fontsize=11)
ax.set_title("(a) Group-specific regression lines", fontsize=12)
ax.legend(fontsize=9); ax.grid(True, alpha=0.3)
# (b) 縮小推定の確認
ax = axes[1]
ax.scatter(u_ols, u_hat, s=n_per_group*3, alpha=0.6, c="steelblue",
edgecolors="black", linewidth=0.5)
lim = max(abs(u_ols).max(), abs(u_hat).max()) + 1
ax.plot([-lim, lim], [-lim, lim], "r--", linewidth=1, label="y = x (no shrinkage)")
ax.plot([-lim, lim], [0, 0], "gray", linewidth=0.5)
ax.plot([0, 0], [-lim, lim], "gray", linewidth=0.5)
ax.set_xlabel("OLS estimate of $u_j$", fontsize=11)
ax.set_ylabel("BLUP of $u_j$", fontsize=11)
ax.set_title("(b) Shrinkage effect (size = $n_j$)", fontsize=12)
ax.legend(fontsize=9); ax.grid(True, alpha=0.3)
# (c) ランダム効果の分布
ax = axes[2]
ax.hist(u_hat, bins=15, density=True, alpha=0.6, color="steelblue",
edgecolor="black", label="BLUP estimates")
x_range = np.linspace(-6, 6, 200)
ax.plot(x_range, 1/(np.sqrt(2*np.pi*sigma_u2_hat))*np.exp(-x_range**2/(2*sigma_u2_hat)),
"r-", linewidth=2, label=f"$N(0, {sigma_u2_hat:.2f})$")
ax.set_xlabel("Random effect $u_j$", fontsize=11)
ax.set_ylabel("Density", fontsize=11)
ax.set_title("(c) Distribution of random effects", fontsize=12)
ax.legend(fontsize=9); ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("mixed_effects.png", dpi=150, bbox_inches="tight")
plt.show()
この結果から、混合効果モデルの重要な特徴が確認できます。
-
グループごとの回帰直線（図a）: 各グループ（学校）が異なる切片を持つ回帰直線が描かれています。黒い太線は全体の固定効果による回帰直線で、各グループの直線はこれをランダム効果 $u_j$ だけ上下にシフトしたものです。傾きは全グループで共通（固定効果）であり、切片のみがグループごとに異なります
-
縮小推定（図b）: BLUPによるランダム効果の推定値がOLS推定値よりもゼロに近い（原点に向かって縮小されている）ことが確認できます。点の大きさはグループサイズ $n_j$ に比例しており、サンプルサイズが小さいグループほど強く縮小されていることがわかります。大きな点（データの多いグループ）は45度線に近く、小さな点（データの少ないグループ）は原点に近い位置にあります
-
ランダム効果の分布（図c）: 推定されたランダム効果のヒストグラムが正規分布の仮定と整合的であることが確認できます。BLUPは縮小推定なので、分布の裾が真の正規分布よりもやや短くなる傾向がありますが、全体的な形状は正規分布に近いです
3つのアプローチの比較
「完全プーリング」（グループ構造を無視）、「プーリングなし」（各グループを独立に推定）、「部分的プーリング」（混合効果モデル）の3つのアプローチを比較します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
n_groups = 20
n_per_group = np.random.randint(3, 25, n_groups)
sigma_u_true = 2.0; sigma_true = 1.5
beta0_true = 5.0
u_true = np.random.randn(n_groups) * sigma_u_true
# データ生成
group_means_true = beta0_true + u_true
y_list = []
group_list = []
for j in range(n_groups):
yj = group_means_true[j] + sigma_true * np.random.randn(n_per_group[j])
y_list.append(yj)
group_list.append(np.full(n_per_group[j], j))
y = np.concatenate(y_list)
group = np.concatenate(group_list).astype(int)
# (1) 完全プーリング: グループ無視、全体平均
grand_mean = y.mean()
# (2) プーリングなし: グループごとの平均
group_means_ols = np.array([y[group==j].mean() for j in range(n_groups)])
# (3) 部分的プーリング: 混合効果モデルのBLUP
sigma2_hat = np.mean([np.var(y[group==j], ddof=1) if n_per_group[j]>1 else sigma_true**2
for j in range(n_groups)])
between_var = np.var(group_means_ols, ddof=1) - sigma2_hat / np.mean(n_per_group)
sigma_u2_hat = max(between_var, 0.01)
beta0_hat = np.sum([n_per_group[j]*y[group==j].mean()/(sigma_u2_hat + sigma2_hat/n_per_group[j])
for j in range(n_groups)]) / \
np.sum([n_per_group[j]/(sigma_u2_hat + sigma2_hat/n_per_group[j])
for j in range(n_groups)])
group_means_blup = np.zeros(n_groups)
for j in range(n_groups):
shrink = sigma_u2_hat / (sigma_u2_hat + sigma2_hat / n_per_group[j])
group_means_blup[j] = beta0_hat + shrink * (group_means_ols[j] - beta0_hat)
# 可視化
fig, ax = plt.subplots(figsize=(10, 6))
idx_sorted = np.argsort(n_per_group)
for rank, j in enumerate(idx_sorted):
ax.plot([rank, rank, rank], [group_means_ols[j], group_means_blup[j], grand_mean],
"k-", linewidth=0.3, alpha=0.3)
ax.scatter(rank, group_means_ols[j], s=n_per_group[j]*5, c="coral",
edgecolors="black", linewidth=0.5, zorder=3)
ax.scatter(rank, group_means_blup[j], s=n_per_group[j]*5, c="steelblue",
edgecolors="black", linewidth=0.5, zorder=3)
ax.scatter(rank, group_means_true[j], s=30, c="black", marker="x", zorder=4)
ax.axhline(grand_mean, color="green", linewidth=1.5, linestyle="--", label="Complete pooling")
ax.scatter([], [], s=50, c="coral", edgecolors="black", label="No pooling (OLS)")
ax.scatter([], [], s=50, c="steelblue", edgecolors="black", label="Partial pooling (BLUP)")
ax.scatter([], [], s=30, c="black", marker="x", label="True group mean")
ax.set_xlabel("Group (sorted by $n_j$)", fontsize=11)
ax.set_ylabel("Group mean", fontsize=11)
ax.set_title("Complete vs No pooling vs Partial pooling", fontsize=12)
ax.legend(fontsize=9, loc="upper left")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("pooling_comparison.png", dpi=150, bbox_inches="tight")
plt.show()
3つのアプローチの比較から、混合効果モデル（部分的プーリング）の利点が明確に読み取れます。
-
完全プーリング（緑の破線）: グループ構造を完全に無視し、全データの平均を使います。すべてのグループに同じ値を割り当てるため、グループ間の真の差異が失われます
-
プーリングなし（赤い点）: 各グループを独立に推定します。サンプルサイズが大きいグループ（右側）は真の値（黒いx）に近い推定を与えますが、サンプルサイズが小さいグループ（左側）は大きくばらつき、極端な値をとることがあります
-
部分的プーリング（青い点）: 混合効果モデルのBLUPは、サンプルサイズの大きいグループではOLS推定値に近く、小さいグループでは全体平均に近い値を取ります。これは最適な妥協であり、データが少ないグループの推定を安定化させつつ、データが多いグループの固有情報を保持します
まとめ
本記事では、混合効果モデルの理論を解説しました。
- モデル: $\bm{y} = \bm{X}\bm{\beta} + \bm{Z}\bm{u} + \bm{\epsilon}$。固定効果 $\bm{\beta}$（未知の定数）とランダム効果 $\bm{u} \sim N(\bm{0}, \bm{G})$（確率変数）を含む
- 周辺分布: $\bm{y} \sim N(\bm{X}\bm{\beta}, \bm{Z}\bm{G}\bm{Z}^\top + \bm{R})$。ランダム効果がグループ内相関を導入する
- ICC: $\sigma_u^2/(\sigma_u^2 + \sigma^2)$ でグループ内の類似性を定量化。0.05以上なら混合効果モデルの使用を推奨
- REML: 分散成分の過小推定を補正する推定法。固定効果の自由度損失を考慮
- BLUP: ランダム効果の最良線形不偏予測量。データの少ないグループは全体平均に縮小される（部分的プーリング）
- 縮小推定: 縮小因子 $\gamma_j = \sigma_u^2/(\sigma_u^2 + \sigma^2/n_j)$ がデータの量に応じて情報源を重み付けする
混合効果モデルは、階層的なデータ構造を持つ多くの実問題に対して、固定効果モデルや完全なランダム効果モデルよりも適切な分析を提供します。
次のステップとして、以下の記事も参考にしてください。