17. 一般化推定方程式（GEE）：クラスター相関データの周辺モデル

GEEの動機：GLMの独立性仮定とクラスターデータの問題

一般化線形モデルは、各観測が互いに独立であるという前提のもとでスコア方程式を構成します。独立性仮定が成立する標準的な状況では、パラメータ $\beta$ の推定量は次のスコア方程式の解として得られます。

$$\sum_{i} \frac{\partial \mu_i}{\partial \beta^T} \phi^{-1} \text{Var}(Y_i)^{-1} (Y_i – \mu_i) = 0$$

しかし繰り返し測定データやクラスターデータでは、この独立性前提が根本的に成立しません。繰り返し測定データとは、同一個体から複数の時点にわたって収集された観測値の集合を指します。クラスターデータは、学校内の生徒や病院内の患者のように、グループ構造をもつ観測単位の集合を指します。両者はクラスター内の観測が互いに相関するという構造的特性を共有します。

独立性を誤って仮定した場合、クラスター内相関が正のときに標準誤差の推定量は真の分散より小さくなる方向に偏り、信頼区間の過小な幅と第1種過誤確率の膨張をもたらします。例えば、同一患者から反復測定された血圧値を独立な観測として扱うと、治療効果の標準誤差が実際より小さく推定され、有意でない効果が有意と判定されるリスクが高まります。

この問題に対し、1986年にLiangとZegerは周辺モデルとしての推定方程式の枠組みを提案しました。周辺モデルは、クラスター内の条件付き分布を完全に特定することなく、集団全体に対する平均的な効果 $E[Y_{ij} | X_{ij}]$ のみを推定対象とします。この枠組みが一般化推定方程式であり、GEEと略されます。個体固有の変動よりも集団全体への介入効果の評価を目的とする分析において基本的な手法として位置づけられます。

GEEの理論的定式化

GEEは、クラスター $i$（$i = 1, \ldots, N$）の観測ベクトル $Y_i = (Y_{i1}, \ldots, Y_{in_i})^T$、周辺平均ベクトル $\mu_i = (\mu_{i1}, \ldots, \mu_{in_i})^T$、および共変量行列 $X_i$ を基本単位とします。各観測の周辺平均は次式で定義されます。

$$E[Y_{ij} | X_{ij}] = \mu_{ij} = g^{-1}(X_{ij}^T \beta)$$

ここで $g$ はリンク関数であり、$\beta$ は集団平均的な効果を表す $p$ 次元の回帰係数ベクトルです。この定式化が周辺モデルの核心であり、ランダム効果や個体固有の変動を明示的にモデル化しない点で条件付きモデルと根本的に異なります。

パラメータの推定はクラスターを単位とした推定方程式を解くことで実行されます。GEEスコア方程式は次の行列形式で表されます。

$$\sum_{i=1}^{N} D_i^T V_i^{-1} (Y_i – \mu_i) = 0$$

ここで $D_i = \partial \mu_i / \partial \beta^T$ は $n_i \times p$ の微分行列であり、平均構造のパラメータに対する感度を表します。作業分散行列 $V_i$ は、周辺分散を対角成分にもつ行列 $A_i$ と作業相関行列 $R(\alpha)$ を用いて定義されます。

$$V_i = A_i^{1/2} R(\alpha) A_i^{1/2}$$

$R(\alpha)$ は推定すべき相関パラメータ $\alpha$ をもつ $n_i \times n_i$ の対称行列であり、クラスター内の観測間相関の近似構造を表します。推定はIRLS（反復重み付き最小二乗法）に類似した反復アルゴリズムで実行されます。$\beta$ の更新ステップと $\alpha$ の更新ステップを交互に繰り返すことで推定値が収束します。

GEEの重要な統計的性質として、作業相関構造 $R(\alpha)$ が真の相関と異なる形式で指定されていても、$\beta$ の推定量は一致性を保持します。これはGEEが準尤度の枠組みと等価であり、平均構造 $E[Y_{ij} | X_{ij}] = \mu_{ij}$ の正しい定式化のみが一致性に要求されることによります。ただし、平均モデルが誤って指定された場合には一致性は保証されません。また、欠測データが後述の完全ランダム欠測の条件を満たすことが、一致推定の別途要件となります。

作業相関構造の種類と選択

作業相関行列 $R(\alpha)$ の構造は分析者が事前に指定します。主要な4種類の構造を以下に示します。

独立構造は $R = I$（単位行列）であり、クラスター内のすべての観測間相関を0と仮定します。推定パラメータを要せず計算負荷が最小であり、後述のSandwich推定量と組み合わせることでロバストな推論が実現されます。交換可能構造は、クラスター内任意の2観測間の相関を一様に $\alpha$ と仮定し、行列の全非対角成分が $\alpha$（$0 \leq \alpha < 1$）となります。推定パラメータは1つであり、対象間の測定順序が定義されない集落研究に適合しやすいです。AR(1)構造は、時間ラグ $|t - s|$ に対して相関が $\rho^{|t-s|}$ と指数的に減衰することを仮定します。縦断研究で時系列的な構造が明確な場合に適用されます。非構造化構造は、クラスター内の全ペアに独立したパラメータを割り当てます。$n_i \times (n_i - 1)/2$ 個の相関パラメータを推定するため柔軟性は最大ですが、パラメータ数も最大となります。

（Fig1. 4種類の作業相関構造（独立・交換可能・AR(1)・非構造化）の相関行列ヒートマップ比較（α = ρ = 0.5、次元5×5））

相関構造の選択にはQIC（準尤度情報量規準）が用いられます。QICは独立作業相関下での準尤度に複雑さの罰則項を加えた指標であり、AICに相当する役割をクラスター相関データの文脈で果たします。候補となる相関構造のなかでQICが最小となる構造を選択することが推奨されます。ただしQICは平均構造の比較には使用せず、相関構造の選択のみを目的とする点でAICとは運用が異なります。

AR(1)構造を適用する場合、各クラスター内での観測の時間順序が明確に定義されていることが前提となります。非構造化構造を採用する場合には、クラスターサイズの均一性が数値的安定性の要件となります。非構造化構造は小標本や不均一なクラスターサイズの状況で数値的に不安定になる点が限界として知られています。作業相関の誤指定は $\beta$ の推定効率を低下させますが、一致性は保持されます。

Sandwich推定量とロバスト標準誤差

GEEにおける標準誤差の推定は、作業相関構造の誤指定に対してロバストなSandwich型の分散推定量が標準的に用いられます。Liang-Zeger型Sandwich分散推定量は次の形式をとります。

$$V_S = A^{-1} B A^{-1}$$

ここで $A = \sum_{i=1}^{N} D_i^T V_i^{-1} D_i$ は期待情報行列に相当し、$B = \sum_{i=1}^{N} D_i^T V_i^{-1} \widehat{\Sigma}_i V_i^{-1} D_i$ は各クラスターの残差外積 $\widehat{\Sigma}_i = (Y_i – \hat{\mu}_i)(Y_i – \hat{\mu}_i)^T$ を積み上げた経験的な項です。この構造からSandwich推定量はHuber-White型の経験分散推定量とも呼ばれ、残差の実測ばらつきに基づいて分散を補正する直感的な解釈をもちます。

モデルベースの推定量（ナイーブ推定量）は $A^{-1}$ の形式であり、作業相関構造が真の共分散と完全に一致している場合にのみ有効な推定量となります。Sandwich推定量は $A^{-1}$ に経験的情報 $B$ を組み込むことで、相関構造の誤指定があっても漸近的に有効な標準誤差を与えます。Independence作業相関とSandwich推定量を組み合わせる戦略は、相関構造に関する仮定を最小化しながらロバストな推論を実現するデフォルト的な選択肢として広く採用されています。

Sandwich推定量の漸近的有効性はクラスター数 $N$ が十分に大きい（$N \to \infty$）ことを理論的前提とします。クラスター数が少ない場合（目安として30未満）、$B$ の経験的推定が不安定となり、Sandwich推定量が真の分散を過小推定する傾向を示します。この問題に対し、HC3型の小標本補正（CR2/CR3と呼ばれる）が有効です。CR3補正はレバレッジ行列を用いて各クラスターの残差項に調整を加えることで、クラスター数が少ない状況でもカバレッジ率を改善します。小標本では自由度補正付きの推定量を積極的に使用することが推奨されます。

（Fig2. クラスター内相関の強さとSandwich推定量・ナイーブ推定量の標準誤差の比較（500反復シミュレーション））

欠測データへの感度とGEEの仮定検証

欠測の発生メカニズムは、欠測指示変数 $M$（観測されていれば0、欠測であれば1）と観測値・欠測値・共変量の条件付き独立性によって分類されます。GEEの文脈での完全ランダム欠測（MCAR）は次式で定義されます。

$$P(M | Y_{obs}, Y_{mis}, X) = P(M | X)$$

この条件は、欠測の発生確率が完全に観測された共変量 $X$ に依存することを許容しますが、観測済みアウトカム $Y_{obs}$ や欠測アウトカム $Y_{mis}$ には依存しないことを要求します。すなわちGEE標準のMCARは、共変量によって欠測が説明される場合には一致性を保持しますが、アウトカム値の大小が欠測確率に影響する場合には一致性が崩れます。MCARの条件のもとでGEEが一致推定量を与える根拠は、推定方程式の各項の期待値がゼロとなることにより不偏性が漸近的に成立するためです。

ランダム欠測（MAR）は $P(M | Y_{obs}, Y_{mis}, X) = P(M | Y_{obs}, X)$ と定義され、欠測が観測済みアウトカム $Y_{obs}$ には依存できるが欠測アウトカム $Y_{mis}$ には依存しない条件を指します。標準GEEはMAR下では残差の期待値がゼロとならず推定方程式の均衡が崩れるため、治療効果の推定に系統的なバイアスが生じます。

この問題への対処として、IPW-GEE（逆確率重み付きGEE）は観測確率の推定モデルを別途構築し、各観測に観測確率の逆数 $\hat{\pi}_{ij}^{-1}$ を重みとして付与することで推定方程式を修正します。IPW-GEEはMAR条件下での一致性を回復する枠組みを提供しますが、欠測非ランダム（MNAR）、すなわち欠測が欠測アウトカム $Y_{mis}$ 自体に依存する場合には、標準GEEおよびIPW-GEEのいずれも一致推定を保証しません。

MCARの検定にはLittle’s testが利用されますが、この検定は一般に検出力が低く、帰無仮説（MCAR）が棄却されなくても実際にはMAR以上の依存構造が存在する可能性を排除できません。そのため、欠測メカニズムに関する仮定に対して感度分析を実施し、推定結果の頑健性を評価することが不可欠です。MAR・MNARデータでの標準GEEの適用は治療効果の一致推定を保証しない点を、結果の報告において明示することが求められます。

臨床試験・縦断研究への応用

GEEは繰り返し測定アウトカムをもつ無作為化比較試験の主解析手法として広く採用されます。2型糖尿病患者を対象とした比較試験において、HbA1c値が複数時点で測定される場合、降糖治療の集団平均的効果は次の線形予測子で定式化されます。

$$\eta_{ij} = \beta_0 + \beta_1 \text{Time}_j + \beta_2 \text{Treat}_i + \beta_3 (\text{Time}_j \times \text{Treat}_i)$$

ここで $\beta_3$ は時間と治療の交互作用効果を表し、治療群と対照群の時間的変化パターンの差を推定します。作業相関構造に交換可能構造を指定しSandwich標準誤差を用いることで、測定間相関に対してロバストな推論が実現されます。2値応答にlogitリンクを用いた場合、GEEの回帰係数の指数変換 $\exp(\hat{\beta}_k)$ は集団平均的なオッズ比を表します。これは特定個体のオッズ比ではなく、対象集団全体における共変量1単位増加に伴う平均的なオッズの変化として解釈されます。集団平均的治療効果はランダム化によって選択バイアスが取り除かれた集団への介入の影響として、政策立案や集団介入の根拠となる推定量を与えます。

クラスター単位データへの適用もGEEの重要な応用領域です。複数の病院を通じた多施設臨床研究では、患者が病院内にネストされたクラスター構造を形成します。このとき各病院の患者群を1クラスターとして、病院間の異質性を交換可能作業相関構造でモデル化することで、病院内相関を適切に調整しながら集団平均的な治療効果を推定できます。この場合のGEE係数は、施設の文脈を平均化した集団全体への介入効果の推定値として解釈されます。

GEEは個体固有の効果を推定する構造を持たないため、個々の患者の治療反応予測や被験者特有の変化軌跡の記述には適していません。また、脱落患者が病状の悪化と関連する場合には欠測がMNARの性質を帯び、標準GEEは治療効果の一致推定量を与えなくなります。この状況ではIPW-GEEと感度分析の併用が不可欠です。周辺効果の推定は集団への政策的介入の評価には有用ですが、個別の臨床意思決定における患者固有のリスク予測には直接対応しない点を解釈に反映する必要があります。

GEEとGLMMの比較：周辺モデルと条件付きモデルの選択指針

GLMM（一般化線形混合モデル）はランダム効果 $u_i \sim N(0, \sigma^2)$ を明示的に組み込み、条件付き平均 $E[Y_{ij} | X_{ij}, u_i] = g^{-1}(X_{ij}^T \beta^{(SS)} + u_i)$ をモデル化します。リンク関数が非線形である場合、GLMMの条件付き平均を周辺化した周辺平均はGEEが直接モデル化する $E[Y_{ij} | X_{ij}]$ と数値的に一致しません。

logistic回帰では、GLMMの条件付き係数 $\beta^{(SS)}$ は対応するGEEの集団平均的係数 $\beta^{(PA)}$ より絶対値で大きくなる傾向を示します。この乖離はランダム効果の分散 $\sigma^2$ が増大するほど拡大します。GLMMの係数は「ランダム効果を固定した条件での個体内変化」を表すのに対し、GEEの係数は「集団全体での平均的な共変量効果」を表します。これらは同一の解釈や数値を共有しません。連続アウトカムにidentityリンクを用いる線形混合モデルでは、GLMMとGEEは数値的に一致する結果を与えます。

（Fig3. GEE（周辺効果）とGLMM（条件付き効果）の係数推定量の比較：ランダム効果分散の増大に伴う乖離（logistic回帰））

集団への介入効果を評価する研究設問では GEEが適切であり、個体内変化の軌跡やランダム効果の大きさを推定する研究設問ではGLMMが適切です。研究設問が「集団全体への介入の影響」を問うものであればGEE、「個体ごとの変化の構造」を問うものであればGLMMを選択することが基本的な指針となります。

GEEは個体レベルの予測やランダム効果の推定を行う構造を持たないため、個体追跡や個別予測が目的の研究には不向きです。GLMMはLaplace近似またはGauss-Hermite数値積分によってランダム効果を周辺化しますが、ランダム効果の次元が高次元になると近似精度が低下するという問題が生じます。