23. 線形混合効果モデル：階層データの回帰

階層データとOLSの限界

通常の最小二乗法は、誤差項が互いに独立かつ等分散であるという仮定のもとで成立します。しかし、患者が病院に属する、生徒が学校に属する、あるいは同一被験者から複数時点で測定が行われる繰り返し測定データでは、同一グループ内の観測値が類似した傾向を示し、この独立性仮定が成立しません。

全分散$\sigma^2_{total}$は、グループ間変動$\sigma^2_b$とグループ内変動$\sigma^2_w$に分解できます。

$$\sigma^2_{total} = \sigma^2_b + \sigma^2_w$$

ICC（クラス内相関係数）は、全変動に占めるグループ間変動の割合として次のように定義されます。

$$\rho = \frac{\sigma^2_b}{\sigma^2_b + \sigma^2_w}$$

$\rho$が0に近い場合はグループ構造の影響が無視できますが、$\rho$が大きい場合は同一グループ内の観測値間の相関が強く、有効なサンプルサイズが総観測数$n$ではなくグループ数$m$に近づきます。たとえば、$\rho = 0.5$で各グループのサイズが10の場合、有効サンプルサイズは観測総数の約17%程度まで低下することが知られています。ICCが高い状況では、最小二乗法が誤差項の相関構造を無視するため、固定効果の標準誤差が過小に推定されます。その結果、第1種過誤率（真の帰無仮説を棄却する確率）が名目水準（通常5%）を超えて膨張し、統計的推測の信頼性が損なわれます。階層データの典型例として、患者（被験者）・病院（施設）・臨床試験の三層構造や、生徒・クラス・学校の二層構造が挙げられます。

線形混合効果モデルの数理的定義

線形混合効果モデル（以下LMM）は、固定効果と変量効果を明示的に組み込むことで、グループ構造を持つデータの分析に適した枠組みを提供します。行列形式での定義は次の通りです。

$$\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{Z}\mathbf{u} + \boldsymbol{\varepsilon}$$

$\mathbf{y}$は$n \times 1$の応答ベクトル、$\mathbf{X}$は$n \times p$の固定効果設計行列、$\boldsymbol{\beta}$は$p \times 1$の固定効果パラメータベクトルです。$\mathbf{Z}$は$n \times q$の変量効果設計行列、$\mathbf{u}$は$q \times 1$の変量効果ベクトル、$\boldsymbol{\varepsilon}$は$n \times 1$の残差ベクトルです。固定効果はすべての観測値に共通する平均的な関係を表し、変量効果はグループごとの個別の偏差を表します。

変量効果と残差の分布はそれぞれ次のように仮定されます。

$$\mathbf{u} \sim \mathcal{N}(\mathbf{0},\, \mathbf{G}), \quad \boldsymbol{\varepsilon} \sim \mathcal{N}(\mathbf{0},\, \mathbf{R})$$

$\mathbf{G}$は変量効果の$q \times q$分散共分散行列、$\mathbf{R}$はグループ内残差の$n \times n$分散共分散行列です。$\mathbf{u}$と$\boldsymbol{\varepsilon}$は独立であると仮定されます。応答ベクトルの周辺共分散行列$\mathbf{V}$は次の式で導出されます。

$$\mathbf{V} = \mathbf{Z}\mathbf{G}\mathbf{Z}^{\top} + \mathbf{R}$$

$\mathbf{y}$の周辺分布はこの共分散構造のもとで次のようになります。

$$\mathbf{y} \sim \mathcal{N}(\mathbf{X}\boldsymbol{\beta},\, \mathbf{V})$$

固定効果$\boldsymbol{\beta}$の推定量はBLUE（最良線形不偏推定量）であり、不偏推定量の中で最小分散を達成します。変量効果$\mathbf{u}$の予測量はBLUP（最良線形不偏予測量）であり、$\mathbf{y}$が観測された条件のもとでの$\mathbf{u}$の条件付き期待値として導出されます。一般化最小二乗法（GLS）は$\mathbf{V}$を既知として固定効果を最適線形不偏推定しますが、LMMはさらに$\mathbf{V}$を変量効果の分散成分から導出し、その分散成分を同時推定する点でGLSを一般化しています。

主要な仮定として、変量効果$\mathbf{u}$の正規性、グループ内残差$\boldsymbol{\varepsilon}$の正規性と独立性、および変量効果と残差ベクトルの独立性が要求されます。

ML推定とREML推定

LMMのパラメータ（固定効果$\boldsymbol{\beta}$および分散成分$\mathbf{G}, \mathbf{R}$）は、最尤法（ML）または制限最尤法（REML）によって推定されます。MLの対数尤度は次の形で表されます。

$$\ell(\boldsymbol{\beta},\mathbf{G},\mathbf{R}) = -\frac{1}{2}\Bigl[n\log(2\pi) + \log|\mathbf{V}| + (\mathbf{y} – \mathbf{X}\boldsymbol{\beta})^{\top}\mathbf{V}^{-1}(\mathbf{y} – \mathbf{X}\boldsymbol{\beta})\Bigr]$$

ML推定では固定効果$\boldsymbol{\beta}$と分散成分を同時推定するため、固定効果の推定に消費した$p$個の自由度が分散成分の推定において考慮されません。その結果、分散成分が過小推定されます（$n-p$ではなく$n$で除することに等価なバイアスが生じます）。グループ数が少ない場合や固定効果の次元が大きい場合に、このバイアスは特に顕著になります。

REMLは、固定効果に直交する$n-p$次元の変換空間に制限したプロファイル対数尤度を最大化することで、この問題を解消します。REML対数尤度は次のように表されます。

$$\ell_{REML}(\mathbf{G},\mathbf{R}) = -\frac{1}{2}\Bigl[\log|\mathbf{V}| + \log|\mathbf{X}^{\top}\mathbf{V}^{-1}\mathbf{X}| + \mathbf{y}^{\top}\mathbf{P}\mathbf{y}\Bigr]$$

ここで$\mathbf{P} = \mathbf{V}^{-1} – \mathbf{V}^{-1}\mathbf{X}(\mathbf{X}^{\top}\mathbf{V}^{-1}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{V}^{-1}$は射影行列です。REMLは分散成分の不偏推定量を与えます。反復アルゴリズムとして、EM法（期待値最大化法）、Newton-Raphson法、Fisher Scoringアルゴリズムが用いられ、収束基準を満たすまで分散成分を逐次更新します。

使い分けの基準として、固定効果構造の比較（例：交互作用項の有無の尤度比検定）にはMLを使用し、分散成分の推定と報告にはREMLを使用します。REMLでは固定効果構造が異なるモデル間の尤度比検定は不適切であり、REMLの対数尤度は異なる設計行列$\mathbf{X}$のもとで比較不可能です。また、グループ数が少ない場合には反復アルゴリズムの収束が不安定になるリスクがあります。

（Fig3. ML vs REML による分散成分推定：グループ数別バイアスと不確実性）

変量効果の構造設計

変量効果の構造はデータの階層性と研究仮説に応じて設計します。ランダム切片のみのモデルでは、グループ$j$の変量効果$u_j$はスカラーであり、$\mathbf{G}$は単一の分散パラメータとなります。

$$\mathbf{G} = \sigma^2_u$$

このモデルはグループごとの切片変動を許容しますが、応答変数と説明変数の関係（傾き）はすべてのグループで共通と仮定します。ランダム切片とランダム傾きの両方を含むモデルでは、変量効果ベクトル$\mathbf{u}_j = (u_{0j},\ u_{1j})^{\top}$が各グループに対して定義され、$\mathbf{G}$は$2 \times 2$の分散共分散行列となります。

$$\mathbf{G} = \begin{pmatrix} \sigma^2_{u_0} & \sigma_{u_0 u_1} \\ \sigma_{u_0 u_1} & \sigma^2_{u_1} \end{pmatrix}$$

$\sigma^2_{u_0}$はグループ間の切片分散、$\sigma^2_{u_1}$は傾き分散、$\sigma_{u_0 u_1}$は切片偏差と傾き偏差の共分散です。共分散が正の場合、切片が高いグループは傾きも大きい傾向があることを意味します。

グループ構造の種類としてネスト効果とクロス効果が区別されます。ネスト効果では一方のグループが他方に完全に包含され（例：患者が病院にネストされる）、クロス効果では両グループの組み合わせが多数の場面で観測されます（例：複数施設で同一の治療群と対照群が設定される場合）。それぞれ異なる設計行列$\mathbf{Z}$として実装されます。

LMMにおける部分プーリングは、ノープーリング（グループごとに独立推定）と完全プーリング（グループを無視した一括推定）の中間に位置します。グループサイズが小さいほど、グループ固有の推定値は全体分布の平均方向に縮小されます。この縮小推定は、データが少ないグループに過大な変動を帰属させることを防ぐための統計的調整機構であり、グループ間で情報を借用することで推定の安定性を高めます。

（Fig1. 完全プーリング・部分プーリング（LMM）・ノープーリングによる回帰直線の比較）

（Fig2. ランダム切片の縮小効果：ノープーリング推定値から全体平均への収縮）

グループ数が少ない場合（特に5グループ未満）は変量効果の分散成分推定が不安定になり、信頼区間の幅が著しく拡大します。また、必要以上にランダム傾きを追加すること（過剰パラメータ化）はモデルの収束失敗やsingular fitの原因となるため、情報量規準を用いた変量効果構造の選択が推奨されます。

モデルの仮定と診断手法

LMMの推測の妥当性は以下の主要仮定が満たされているかどうかに依存します。変量効果$\mathbf{u}$の正規性、グループ内残差$\boldsymbol{\varepsilon}$の独立性と等分散性、そして変量効果と固定効果説明変数の直交性（内生性の非存在）の3点が特に重要です。

変量効果の正規性はQQプロットによって診断されます。推定されたBLUP値を正規スコアに対してプロットし、直線からの顕著な逸脱は正規性仮定の違反を示唆します。グループ内残差については、条件付き残差（フィット値と残差のプロット）でランダムな散らばりを確認します。残差のフィット値依存性や周期的パターンは、残差の等分散性または系列独立性の違反を示します。

変量効果と固定効果説明変数の内生性が存在する場合（例：グループ属性と応答変数の共通未測定要因が存在する場合）、固定効果の推定量にバイアスが生じます。この問題はハウスマン検定などで検出できます。影響観測値の評価には、OLSのCook距離をLMMに拡張した指標が用いられ、特定観測値の除外がBLUEおよびBLUPの推定値に与える変化量として定義されます。

共分散構造の誤設定（例：残差に自己相関が存在するのに$\mathbf{R} = \sigma^2\mathbf{I}$を仮定する場合）は、固定効果の標準誤差を歪め、推測の信頼性を損ないます。変量効果の分布が非正規の場合（例：二峰性や歪みが存在する場合）、REMLの分散成分推定は完全な頑健性を持たず、特にサンプルサイズが小さい状況では推測精度が低下します。共分散構造の候補として、複合対称性、AR(1)、非構造化行列などをAIC（赤池情報量規準）やBIC（ベイズ情報量規準）で比較することが診断の一助となります。残差のQQプロットとスケール-ロケーションプロットを組み合わせた診断により、分散均一性と正規性の両面を同時に確認することが推奨されます。

生物統計学への応用：臨床試験の反復測定データ

臨床試験において、複数の時点で測定が繰り返される縦断データは自然な階層構造を持ちます。被験者$i$の時点$t$での応答変数$y_{it}$に対して、被験者をグループ単位とするランダム切片モデルを当てはめることで、治療効果と被験者間変動を分離して推定できます。固定効果として治療群ダミー変数、時点変数、および時点×治療の交互作用項を含めることで、治療効果の時間的推移を推定できます。時点×治療の交互作用係数は、治療群間で変化率が統計的に異なるかどうかを評価する固定効果パラメータとして解釈されます。

LMMが繰り返し測定ANOVAより優れる主要な点として、欠測データ処理の枠組みが挙げられます。LMMは直接最大尤度法を利用するため、欠測がMAR（ランダム欠測）の仮定のもとで発生していれば、欠測被験者を自動的に除外する完全ケース解析と比較して偏りの少ない固定効果推定が可能です。繰り返し測定ANOVAは全時点のデータが揃った被験者のみを解析対象とするため、MARが成立する場面でも選択バイアスが生じます。また、時点間の分散均一性（球面性仮定）が不要な点も、LMMの実用的な利点です。

治療効果は固定効果$\boldsymbol{\beta}$中の治療群係数として推定され、被験者間変動はランダム切片の分散$\sigma^2_{u_0}$として定量化されます。GLSとの比較では、LMMはグループ（被験者）レベルの変量効果を明示的に推定するため、個人差の分量と治療効果を明確に分離できる点が優れています。

適用における主要な仮定と限界として、MAR仮定が成立しないMNAR（非ランダム欠測）型の脱落（例：治療効果が不良な被験者が試験を中断する場合）には対応できません。この場合、感度分析や共有パラメータモデルの導入が必要です。被験者数や時点数が多い大規模試験では、$\mathbf{V}$の逆行列演算に伴う計算コストが増大します。また、被験者数が少ない小規模試験では分散成分推定が不安定になるため、点推定値だけでなく推定の不確実性を評価することが重要です。

比較モデルと拡張：GLMMと階層ベイズ

グループ構造を持つデータに対するモデル選択の枠組みとして、完全プーリング・部分プーリング・ノープーリングの三分類が有効です。完全プーリングはグループ構造を無視してOLSで一括推定する方法であり、グループ間変動を説明できない代わりに推定の分散が小さいです。ノープーリングはグループごとに独立した固定効果モデルを推定する方法であり、グループ固有の効果を捕捉できますが、グループサイズが小さいと推定分散が大きくなりグループ間の情報が活用できません。LMMによる部分プーリングはこの中間に位置し、グループ間変動の分散構造を利用して各グループの推定値を安定化させます。

応答変数が二値（ロジスティック回帰に相当）やカウントデータ（ポアソン回帰に相当）の場合、LMMを一般化線形混合効果モデル（GLMM）に拡張できます。GLMMでは線形予測子$\mathbf{X}\boldsymbol{\beta} + \mathbf{Z}\mathbf{u}$にリンク関数を通じて非正規応答を扱います。GLMMの周辺尤度は解析的に求まらないため、Laplace近似またはGauss-Hermite求積などの数値積分近似を必要とし、近似精度が推定結果に影響します。

LMMと階層ベイズモデルの関係として、変量効果に対して弱情報事前分布（$\mathbf{u} \sim \mathcal{N}(\mathbf{0}, \mathbf{G})$に対応する半正規事前分布など）を設定した場合、階層ベイズモデルはLMMと実質的に等価な推定を与えます。ベイズ枠組みでは、MCMC法（マルコフ連鎖モンテカルロ法）による事後分布からの直接サンプリングが可能なため、小標本での不確実性評価において優れた柔軟性を持ちます。個体成長曲線や薬物動態モデルのように、応答変数が時間に対して非線形な関係を持つ場合には、非線形混合効果モデルへの拡張が必要となります。階層の深さが増すほど周辺尤度の計算負荷が指数的に増大するため、深い階層構造を持つモデルでは計算効率の高い近似手法の選択が重要になります。