GEE（一般化推定方程式）とは ― 周辺モデルの考え方とRで学ぶ反復測定データ解析 ―
この記事でわかること
・GEE（一般化推定方程式）とは何か ― GLMを「相関のある反復測定データ」へ拡張する枠組み
・周辺モデル（集団平均）と条件付きモデル（個人内）の違い ― 同じデータでも推定している対象が違う
・作業相関構造（working correlation）と頑健分散（サンドイッチ分散）の考え方と、その頑健性の意味
・Rの
geepack パッケージによる実装 ― 二値アウトカムの反復測定データを実際に解析する・GEEとGLMM・MMRMの使い分け ― 臨床試験のどの場面でどれを選ぶか
はじめに
臨床試験のデータを扱っていると、「1人の被験者から1つの観測値」というきれいな状況はむしろ少数派だと気づきます。ベースライン、4週、8週、12週……と同じ患者を繰り返し測定するのが当たり前ですし、多施設試験では同じ施設の患者同士が似た傾向を持つこともあります。このとき、観測値は独立ではありません。
この相関を無視して通常のロジスティック回帰やポアソン回帰を当てはめると、回帰係数の点推定値はさほど変わらないことが多い一方で、標準誤差が誤ります。直感的には、同じ患者の4回の測定を「独立した4人分の情報」とみなしてしまうため、実際より多くの情報を持っていると錯覚し、信頼区間が狭くなりすぎp値が小さく出すぎる、という方向に歪みます。実務では、有意でないものを有意と誤って主張してしまうリスクに直結します。
連続量のアウトカムを繰り返し測定した場合、製薬業界での定番はMMRM（反復測定混合モデル）で、実装についてはmmrmをRで実装するで詳しく扱いました。しかし、アウトカムが二値（改善した／しない）やカウント（増悪回数など）である場合、正規分布ベースの枠組みはそのままでは使えません。ここで主役になるのがGEE（Generalized Estimating Equations、一般化推定方程式）です。
この記事は、製薬企業の生物統計担当者や、R・SASで反復測定データの解析に取り組む実務家の方を主な読者として想定しています。まずGEEの核心である「周辺モデル」を整理し、続いて推定方程式・作業相関構造・頑健分散という三つの柱を押さえます。その上でRの geepack パッケージで二値アウトカムの反復測定データを解析し、最後にGLMM・MMRMとの使い分けを、同じデータの数値を並べながら確認していきます。
GEEとは ― 「集団平均」を推定する周辺モデル
GEEは、LiangとZegerが1986年に提案した手法です。一言でいえば、GLM（一般化線形モデル）を、相関のある反復測定データやクラスタデータへ拡張したものにあたります。ロジットリンクを使えば二値アウトカムに、対数リンクを使えばカウントデータに、という具合に、GLMで慣れ親しんだリンク関数と分散関数の枠組みをそのまま持ち込めるのが大きな魅力です。
そして、GEEを理解する上で最も重要な概念が周辺モデル（marginal model）です。周辺モデルが推定するのは「集団全体で平均的にどうか」という効果で、たとえば「実薬群の状態良好割合は、プラセボ群と比べて集団平均としてどれだけ高いか」という問いに答えます。これに対して、GLMM（一般化線形混合モデル）のようなランダム効果モデルは条件付きモデル（conditional model）であり、「同じ個人の中で、その人が実薬を受けたときと受けなかったときでどれだけ違うか」という、個体を固定した上での効果を推定します。GLMMそのものについてはGLMM（一般化線形混合モデル）とはで解説しています。
この違いは、単なる言葉づかいの問題ではありません。臨床試験の文脈では「この治療を集団に適用したとき、平均的にどれだけ改善が見込めるか」という問いに答えたい場面が非常に多く、そこでは周辺解釈のほうが自然でストレートです。ICH E9(R1)のエスティマンドの言葉を借りるなら、何を推定したいのかを先に決め、それに合う手法を選ぶという順序が本来の姿であり、GEEとGLMMの選択もまさにその一例です。
GEEのもうひとつの特徴は、完全な尤度を仮定しない準尤度（quasi-likelihood）的なアプローチである点です。アウトカムの同時分布を正しく書き下す必要はなく、必要なのは平均構造と分散関数、そして相関構造の当たりをつけた形だけ。しかもその相関構造は作業（working）相関構造、つまり「とりあえずの作業仮説」として置けばよく、多少外れても致命傷にはなりません。この気楽さこそがGEEの実務的な強みです。
GEEの数理 ― 推定方程式・作業相関構造・頑健分散
推定方程式
名前が示すとおり、この手法の中心にあるのは「尤度の最大化」ではなく「推定方程式を解くこと」です。クラスタ（ここでは被験者）を \(i = 1, \dots, K\) と添字づけると、回帰係数 \(\boldsymbol{\beta}\) は次の方程式の解として定義されます。
\[ \sum_{i=1}^{K} \mathbf{D}_i^{\top} \mathbf{V}_i^{-1} (\mathbf{Y}_i – \boldsymbol{\mu}_i) = \mathbf{0} \]
ここで登場する三つの部品を、日本語でかみくだいておきます。\(\mathbf{Y}_i – \boldsymbol{\mu}_i\) は残差で、被験者 \(i\) について実際に観測されたアウトカムのベクトル \(\mathbf{Y}_i\) と、モデルが予測する平均 \(\boldsymbol{\mu}_i\) とのずれを表します。GEEは、このずれの重み付き合計がちょうどゼロになるように \(\boldsymbol{\beta}\) を選ぶ、という発想です。
\(\mathbf{D}_i\) は平均 \(\boldsymbol{\mu}_i\) を \(\boldsymbol{\beta}\) で偏微分した行列で、「\(\boldsymbol{\beta}\) を少し動かしたら予測平均がどれだけ動くか」という感度を表します。残差の情報を係数の空間へ翻訳する役割だと考えると分かりやすいでしょう。
そして \(\mathbf{V}_i\) が作業共分散行列、すなわち「被験者 \(i\) の4つの観測値は、互いにこのくらい相関しているだろう」という作業仮説です。これが逆行列 \(\mathbf{V}_i^{-1}\) として重みに入ることで、相関の強い観測値の情報を重複して数えない重み付けが自動的にかかります。この \(\mathbf{V}_i\) は
\[ \mathbf{V}_i = \mathbf{A}_i^{1/2} \mathbf{R}_i(\alpha) \mathbf{A}_i^{1/2} \]
と分解されます。\(\mathbf{A}_i\) は各時点の分散（二項分布なら \(\mu(1-\mu)\)）を並べた対角行列、\(\mathbf{R}_i(\alpha)\) が次に述べる作業相関行列です。この分散関数がどこから来るのかは一般化線形モデル(GLM)の“裏側”で指数型分布族の側から整理しています。\(\mathbf{R}_i\) を単位行列にすれば、この式はGLMのスコア方程式そのものに戻ります ― GEEがGLMの素直な拡張だというのはこういう意味です。
作業相関構造
作業相関行列 \(\mathbf{R}_i(\alpha)\) をどう置くかが、GEEを使う際の実質的な設計判断になります。代表的な4つを整理します。
| 構造名 | \(R_i\) の考え方 | 推定する相関パラメータ数 | 向いている場面 |
|---|---|---|---|
| independence（独立） | 単位行列。クラスタ内の相関を一切仮定しない | 0個 | 構造に見当がつかないとき。頑健分散と組み合わせれば最低限の妥当性は確保される |
| exchangeable（交換可能） | どの2観測の相関もすべて同じ \(\alpha\)。順序を入れ替えても同じ | 1個 | クラスタ内に順序がない場合（施設、家族、患者クラスタ）。クラスタランダム化試験の定番 |
| AR(1)（1次自己回帰） | 時点差 \(k\) の相関が \(\alpha^{k}\)。離れるほど相関が減衰する | 1個 | 等間隔で測定される縦断データ。時間が離れるほど関連が弱まる臨床試験の来院データ |
| unstructured（無構造） | すべての時点ペアに別々の相関を割り当てる。仮定を置かない | \(T(T-1)/2\) 個（4時点なら6個） | 時点数が少なく被験者数が多い試験。相関の形を事前に決め打ちしたくないとき |
実務的な指針としては、まず測定の性質から選ぶのが筋のよいやり方です。臨床試験の時点間相関は時間が離れるほど弱くなることが多いので、来院ごとの縦断データならAR(1)が自然な第一候補になります。一方、クラスタ内に順序という概念がない ― 施設、家族、同一診療科の患者群といった場合は交換可能が理にかなっています。時点数が4つ程度と少なく被験者数が十分にある試験なら、無構造を選んで相関の形をデータに語らせるのも現実的です。
頑健分散（サンドイッチ分散推定量）
「作業相関構造の選択を間違えたらどうなるのか」と不安になった方もいるかもしれませんが、実はそこにGEE最大の売りがあります。GEEでは標準誤差を、モデルの仮定に依存する形ではなく頑健分散（robust variance、サンドイッチ分散推定量）として計算します。その性質を一言でいえば ― 作業相関構造が間違っていても、平均構造さえ正しく指定されていれば、回帰係数の推定は一致性を保ち、標準誤差も妥当なものが得られる。これが頑健性（robustness）と呼ばれる性質で、相関構造は「作業仮説」で構わないと先に述べたのはこの裏づけがあるからです。
名前の由来は式の見た目にあります。分散推定量が \(B^{-1} M B^{-1}\) という形をしており、中央の \(M\)（実際の残差から計算される部分、いわば具のハム）を両側から \(B^{-1}\)（モデルの仮定に基づく部分、パン）で挟んでいるように見えるためです。真ん中に実際のデータのばらつきそのものが入っているからこそ、外側のパン（作業相関の仮定）が多少ずれても全体として辻褄が合う、と理解しておくと感覚がつかめます。
頑健分散は万能ではありません。最大の落とし穴は、クラスタ数が少ないときに頑健分散が過小に推定されることです。頑健性は「クラスタ数が十分大きいとき」に成り立つ漸近的な性質であり、クラスタが少ないと中央の \(M\) を推定するための情報が足りず、標準誤差が小さく出て第一種の過誤が膨らみます。目安としてクラスタ数が40程度を下回る場合は、Mancl-DeRouenやKauermann-Carrollといったバイアス補正、あるいは自由度を調整するスモールサンプル補正の検討が必要です。少数施設のクラスタランダム化試験や小規模な第II相試験では、特に注意してください。
RでGEEを実装する ― geepackで反復測定データを解析する
ここからは実データでGEEを動かします。使うのは geepack 同梱の respiratory データで、呼吸器疾患を対象に実施された本物のランダム化比較試験です。2施設・111人（プラセボ群 P が57人、実薬群 A が54人）を各4回の受診で追跡しており、1人4行で計444行という典型的な反復測定の形をしています。アウトカム outcome は各時点で状態が良好なら1をとる二値変数です。二値アウトカムを患者単位で繰り返し測定している——まさにGEEの出番です。
| 変数 | 内容 |
|---|---|
| center / id | 実施施設（1・2）と患者番号。idは施設ごとに1から振り直されている |
| treat | 治療群（P＝プラセボ、A＝実薬） |
| sex / age | 性別（F/M）・年齢。ベースライン共変量 |
| baseline | 試験開始時の呼吸器状態（0＝不良、1＝良好） |
| visit / outcome | 受診時点（1〜4）と、その時点で状態が良好なら1をとるアウトカム |
データの準備とクラスタIDの作成
まずデータを解析できる形に整えます。GEEではこの下準備が結果そのものを左右します。
# geepackの読み込みとデータの準備
library(geepack)
data(respiratory)
resp <- respiratory
# idは施設ごとに1から振り直されているため、施設×患者で一意のクラスタIDを作る
resp$id2 <- interaction(resp$center, resp$id)
# 受診時点は連続量ではなく因子として扱う（時点ごとの効果を自由に推定するため）
resp$visitf <- factor(resp$visit)
# 参照群をプラセボ（P）にする。これで係数treatAが「実薬 vs プラセボ」になる
resp$treat <- relevel(factor(resp$treat), ref = "P")
# GEEは「同じidの行が隣り合っている」ことを前提とするため、クラスタ内で時点順に並べ替える
resp <- resp[order(resp$id2, resp$visit), ]
# 観測された「状態良好」の割合を treat × visit で確認する
round(tapply(resp$outcome, list(resp$treat, resp$visitf), mean), 3)
1 2 3 4
P 0.491 0.386 0.456 0.439
A 0.685 0.704 0.722 0.611
モデルを当てはめる前に生の割合を眺めるのが実務の鉄則です。実薬群Aは 0.685 → 0.704 → 0.722 → 0.611 と推移し、プラセボ群Pの 0.491 → 0.386 → 0.456 → 0.439 を4時点すべてで上回っています。差は20ポイント前後で、時点による逆転もありません。GEEが推定するのはまさにこの群ごとの平均的な割合の差であり、集計表の世界と地続きです。
このデータの
id は施設ごとに1から振り直されています。「施設1の患者3番」と「施設2の患者3番」は別人なのに同じ値です。何も考えずに id = id と指定すると別人の測定が同一患者としてまとめられ、相関の推定も標準誤差も丸ごと狂います。interaction(center, id) で施設×患者の一意IDを作ることが必須です。もうひとつ、geeglm() は同じidの行が隣り合っていることを前提とします。行が散らばっていると作業相関構造（特にAR(1)や無構造）が意図しない順序で当てはめられるため、解析前に必ず order(id, time) でソートしてください。
交換可能相関でGEEを当てはめる
主解析を当てはめます。作業相関構造には交換可能相関（exchangeable、同一患者内のどの2時点の相関も等しいと仮定する構造）を採用します。
# モデル式：治療群・受診時点・ベースライン状態・性別・年齢で調整
f <- outcome ~ treat + visitf + baseline + sex + age
# 交換可能相関を作業相関構造としたGEE（主解析）
fit_exch <- geeglm(f, family = binomial, data = resp,
id = id2, corstr = "exchangeable")
summary(fit_exch)
Call:
geeglm(formula = f, family = binomial, data = resp, id = id2,
corstr = "exchangeable")
Coefficients:
Estimate Std.err Wald Pr(>|W|)
(Intercept) -3.89e-01 6.59e-01 0.35 0.55511
treatA 1.26e+00 3.50e-01 13.01 0.00031 ***
visitf2 -2.40e-01 2.46e-01 0.95 0.32980
visitf3 -1.30e-16 2.64e-01 0.00 1.00000
visitf4 -3.35e-01 2.45e-01 1.86 0.17240
baseline 1.98e+00 3.27e-01 36.60 1.5e-09 ***
sexM -2.77e-01 4.20e-01 0.43 0.51071
age -1.34e-02 1.33e-02 1.01 0.31494
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation structure = exchangeable
Estimated Scale Parameters:
Estimate Std.err
(Intercept) 0.977 0.152
Link = identity
Estimated Correlation Parameters:
Estimate Std.err
alpha 0.364 0.0842
Number of clusters: 111 Maximum cluster size: 4
treatA は係数 1.26、頑健標準誤差 0.350、Wald統計量 13.01、p値 0.00031 で、実薬群が有意に良好という結果です。この
Std.err がすでに頑健標準誤差（サンドイッチ分散に基づく標準誤差）である点に注意してください。geeglm() は既定で頑健分散を返します。alpha＝0.364（標準誤差 0.0842）は同じ患者の異なる受診時点どうしの相関がおよそ0.36あることを意味し、これを無視して444件の独立な観測として扱う危険が数字から見て取れます。
最終行の Number of clusters: 111 Maximum cluster size: 4 も重要です。111人×最大4回という構造どおりにクラスタが認識できています。GEEを回したらまずこの行で検算する癖をつけると事故を減らせます。
オッズ比に変換して解釈する
ロジットリンクなので係数はそのままでは対数オッズ比です。頑健標準誤差を使って、オッズ比と95%信頼区間の形に直します。
# 係数と頑健標準誤差からオッズ比・95%信頼区間を計算する
est <- summary(fit_exch)$coefficients
or <- data.frame(
OR = exp(est[, "Estimate"]),
LCL = exp(est[, "Estimate"] - 1.96 * est[, "Std.err"]),
UCL = exp(est[, "Estimate"] + 1.96 * est[, "Std.err"]),
p = est[, "Pr(>|W|)"]
)
print(or, digits = 3)
OR LCL UCL p
1 0.678 0.186 2.47 0.5551
2 3.531 1.779 7.01 0.0003
3 0.787 0.486 1.27 0.3298
4 1.000 0.596 1.68 1.0000
5 0.716 0.442 1.16 0.1724
6 7.215 3.804 13.69 0.0000
7 0.758 0.333 1.73 0.5107
8 0.987 0.961 1.01 0.3149
行の順番はモデルの係数の並びと同じで、2行目が treatA、6行目が baseline に対応します。
treatA のオッズ比は 3.531（95%信頼区間 1.779–7.01）です。ここで言葉づかいが決定的に重要になります。GEEが推定するのは周辺（marginal）効果ですから、正しい読み方は「実薬群では平均的に、状態が良好となるオッズがプラセボ群の約3.5倍」——すなわち集団平均（population-averaged）としての効果です。「その患者個人のオッズが3.5倍になる」という言い換えは誤りで、この違いは後半のGLMMとの比較ではっきり現れます。
baseline のオッズ比は 7.215（3.804–13.69）と圧倒的に大きく、試験開始時の状態がその後を最も強く予測する因子だとわかります。ベースライン調整で精度が上がるのは臨床試験共通の原則で、その考え方と規制上の位置づけは共変量調整（ANCOVA）徹底解説で整理しています。
時点効果はいずれも有意ではありません。なかでも visitf3 の係数が -1.30e-16、オッズ比がちょうど 1.000 なのは偶然ではありません。実際に数えると、状態良好であった人数は visit 1 が65人、visit 3 も65人とまったく同数です（いずれも111人中、割合0.5856）。因子として投入した時点効果は各水準の周辺度数を再現するように推定されるため、参照水準である visit 1 とイベント数が一致すれば、係数は数値誤差の範囲で厳密にゼロになります。4回の受診を通じて状態が系統的に改善・悪化する傾向は見られない、というのがここでの読み取りです。
頑健分散はなぜ必要か ― クラスタリングを無視したglmとの比較
頑健分散のありがたみは、あえて相関を無視した解析と並べると一目瞭然です。同じモデル式を、クラスタ構造を一切考慮しない通常のロジスティック回帰で当てはめてみます。
# クラスタリングを完全に無視した通常のロジスティック回帰
# （444行すべてを独立な観測として扱ってしまう）
fit_glm <- glm(f, family = binomial, data = resp)
summary(fit_glm)$coefficients
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.95e-01 0.47116 -8.38e-01 4.02e-01
treatA 1.29e+00 0.23349 5.52e+00 3.41e-08
visitf2 -2.41e-01 0.31069 -7.75e-01 4.38e-01
visitf3 3.29e-16 0.31180 1.06e-15 1.00e+00
visitf4 -3.36e-01 0.31064 -1.08e+00 2.79e-01
baseline 2.01e+00 0.23480 8.54e+00 1.30e-17
sexM -2.73e-01 0.28694 -9.50e-01 3.42e-01
age -1.38e-02 0.00856 -1.61e+00 1.08e-01
通常のglmの標準誤差と、同じ独立作業相関を指定したGEEの頑健標準誤差を並べてみます。点推定値が完全に一致する独立作業相関どうしで比べることで、標準誤差だけの違いを純粋に取り出せます。
| 共変量 | タイプ | 通常のglmのSE | GEEの頑健SE | 比（GEE / glm） |
|---|---|---|---|---|
| treat | 被験者間で一定 | 0.23349 | 0.35190 | 1.507（拡大） |
| baseline | 被験者間で一定 | 0.23480 | 0.32835 | 1.398（拡大） |
| age | 被験者間で一定 | 0.00856 | 0.01340 | 1.565（拡大） |
| visitf2 | 被験者内で変化 | 0.31069 | 0.24774 | 0.797（縮小） |
| visitf4 | 被験者内で変化 | 0.31064 | 0.24723 | 0.796（縮小） |
まず点推定値がまったく動いていません。glmの treatA も独立作業相関のGEEも 1.28858 で、小数第5位まで完全に一致します。作業相関を独立と置くとGEEの推定方程式は通常のロジスティック回帰のスコア方程式そのものに戻るため、これは偶然ではなく数学的な必然です。変わるのは標準誤差だけ——頑健分散は「推定値を直すもの」ではなく「不確実性の見積もりを直すもの」なのです。
注目すべきは、その変化が一方向ではないことです。被験者間で一定の共変量（treat・baseline）では、glmがSEを小さく見積もりすぎます。治療群は4回とも変わらないので、同じ患者の4行は治療効果について同じ情報を繰り返しているだけです。それをglmは独立な444人分として数え、実際には111人分しかない情報量を水増しして精度を過大評価します。treat が 0.233 → 0.352 と約1.5倍に広がったのは、その補正分です。
逆に、被験者内で変化する共変量（visit）ではglmがSEを大きく見積もりすぎます。時点間の比較は同じ患者の中での比較で、体質や重症度といった個人差が相殺されるため、むしろ精度は上がります。実際 visitf2 は 0.311 → 0.248、visitf4 は 0.311 → 0.247 と狭くなりました。
表全体を眺めると、この対比はいっそう鮮やかです。被験者間で一定の共変量（treat・baseline・age）はそろって比が1.40〜1.57と拡大し、被験者内で変化する共変量（visitf2・visitf4）はそろって0.80前後へ縮小しています。共変量の性質によって補正の向きがきれいに二分されているわけです。
この非対称性を知らないと「クラスタを無視すると常に有意になりやすい」と誤解します。正しくは「クラスタ内で一定の共変量は過小評価、変わる共変量は過大評価」です。臨床試験で主要な関心となる治療群は前者にあたるため、クラスタリングの無視は偽陽性の方向に効きます。頑健分散が規制対応上の必須要件になっているのは、このためです。
作業相関構造の選び方 ― QICと頑健性
作業相関構造をどう選ぶべきかは、GEEを使い始めた誰もが直面する疑問です。代表的な4構造を当てはめて、その選択が結論をどれだけ左右するのかを数字で確かめます。
代表的な4構造を当てはめて比較する
independence（独立）・exchangeable（交換可能）・ar1（1次自己回帰）・unstructured（無構造）を順に当てはめ、treatA の推定値を並べます。
# 4つの作業相関構造を順に当てはめる
structures <- c("independence", "exchangeable", "ar1", "unstructured")
fits <- lapply(structures, function(cs) {
geeglm(f, family = binomial, data = resp, id = id2, corstr = cs)
})
names(fits) <- structures
# treatA（実薬 vs プラセボ）の推定値・頑健標準誤差・オッズ比を並べる
comp <- do.call(rbind, lapply(structures, function(cs) {
co <- summary(fits[[cs]])$coefficients["treatA", ]
data.frame(corstr = cs,
Est = co[["Estimate"]],
SE = co[["Std.err"]],
OR = exp(co[["Estimate"]]))
}))
print(comp, digits = 5)
corstr Est SE OR
1 independence 1.2886 0.35190 3.6276
2 exchangeable 1.2617 0.34978 3.5313
3 ar1 1.1701 0.35247 3.2223
4 unstructured 1.2477 0.34867 3.4824
推定された相関パラメータも見ておきましょう。無構造では時点の全組み合わせについて個別の相関 \( \alpha \) が推定されます。
# 無構造で推定された時点ペアごとの相関パラメータ
summary(fits[["unstructured"]])$corr
Estimate Std.err
alpha.1:2 0.342 0.1111
alpha.1:3 0.219 0.0984
alpha.1:4 0.345 0.1079
alpha.2:3 0.439 0.1300
alpha.2:4 0.417 0.1319
alpha.3:4 0.415 0.1231
4構造のオッズ比は 3.6276・3.5313・3.2223・3.4824 と3.22〜3.63の範囲に収まり、標準誤差も 0.34867〜0.35247 とほぼ同じです。どの構造を選んでも「実薬群が有意に良好」という結論は変わりません——これこそが頑健分散の恩恵で、「作業（working）相関」という呼び名自体が「作業用の仮の仮定でよい」という含意を持っています。
無構造の相関は alpha.1:2＝0.342 から alpha.2:3＝0.439 まで、おおむね0.2〜0.44の範囲にあり、交換可能相関の alpha＝0.364 はその中央付近です。一方、AR(1)の alpha＝0.501 が想定する「離れるほど相関が弱くなる」パターンはきれいには成り立っていません。1:3 が 0.219 と低いのに 2:4 は 0.417 と高く、受診間隔が等間隔で時間的な減衰がはっきりしないデータではAR(1)が窮屈な仮定になり得るという実例です。
QICによる構造の選択
構造の選択規準としては、AICを準尤度（quasi-likelihood）に拡張したQIC（Quasi-likelihood under the Independence model Criterion）が広く使われます。geepack では QIC() で計算できます。
# QICで4つの構造を比較する（小さいほど良い）
print(t(sapply(fits, QIC)), digits = 3)
QIC QICu Quasi Lik CIC params QICC
independence 518 509 -247 12.6 8 520
exchangeable 518 509 -247 12.4 8 520
ar1 519 510 -247 12.5 8 520
unstructured 518 509 -247 12.4 8 522
QICはAICの準尤度版で、小さいほど良いと判断します。independence・exchangeable・unstructured がいずれも518で並び、ar1 が519とわずかに劣るだけでした。有効パラメータ数に相当するCICも 12.4〜12.6 で、差はごくわずかです。
つまりこのデータでは構造選択の影響が小さいということです。判断材料は当てはまりよりも解釈のしやすさになり、パラメータが \( \alpha \) ひとつだけで「同じ患者の測定どうしは一様に約0.36相関している」と一言で説明できる交換可能相関を主解析に据えるのが妥当と結論づけられます。
ただしQICは、今回のように差が1〜2しかないときには決め手になりません。「QICが最小だったから」という説明は根拠として機能しない、と正直に認識しておくべきです。
「構造の選択は結論を左右しない」というのは、データを見てから選んでよいという意味ではありません。臨床試験、とりわけ検証的試験では作業相関構造を解析計画書（SAP）で事前に固定するのが原則です。複数試して最も都合のよい結果を採れば実質的な多重性が生じ、第一種の過誤が膨らみます。SAPに主解析の構造（多くは交換可能相関）を明記し、他は感度解析として事前に規定しておくのが標準的な進め方です。今回のように結果が安定していれば、「結論が相関の仮定に依存しない」という強い支持材料として報告できます。
GEEとGLMM・MMRMの使い分け ― 周辺効果と条件付き効果
ここまでGEEの推定と作業相関構造を見てきましたが、実務でもっとも質問が多いのは「同じ反復測定データなのに、GEEとGLMM、そしてMMRMのどれを使えばよいのか」という点です。この問いに答えるには、それぞれが推定している対象（エスティマンド）が違うという事実を、数値で腹落ちさせるのが近道です。
そこで、これまでとまったく同じ respiratory データ・まったく同じ共変量に、今度はGLMM（一般化線形混合モデル）をあてはめてみます。患者ごとの体質やベースラインの重症度といった「その人固有のかかりやすさ」を、患者ごとのランダム切片 (1|id2) として明示的にモデルに組み込む、という発想です。GLMMそのものの考え方についてはGLMM（一般化線形混合モデル）とはで詳しく解説していますので、あわせてご覧ください。
library(lme4)
# ランダム切片 (1|id2) で患者ごとの「かかりやすさ」を明示的にモデル化する
# 既定のNelder-Mead法では収束警告が出るため、最適化アルゴリズムをbobyqaに変更している
fit_glmm <- glmer(outcome ~ treat + visitf + baseline + sex + age + (1|id2),
data = resp, family = binomial,
control = glmerControl(optimizer = "bobyqa"))
summary(fit_glmm)
Random effects:
Groups Name Variance Std.Dev.
id2 (Intercept) 3.96 1.99
Number of obs: 444, groups: id2, 111
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -8.04e-01 1.04e+00 -0.77 0.43859
treatA 2.12e+00 5.52e-01 3.84 0.00012 ***
visitf2 -3.84e-01 3.94e-01 -0.98 0.32883
visitf3 4.12e-07 3.94e-01 0.00 1.00000
visitf4 -5.36e-01 3.95e-01 -1.36 0.17399
baseline 3.31e+00 5.96e-01 5.56 2.6e-08 ***
sexM -4.46e-01 6.68e-01 -0.67 0.50461
age -1.99e-02 1.96e-02 -1.01 0.31180
---
まず treatA の係数を見比べてください。GEE（交換可能相関）では 1.26、オッズ比に直すと OR = 3.53 でした。ところが同じデータ・同じ共変量のGLMMでは 2.12、OR = 8.37 です。オッズ比が2倍以上違います。
初めてこの現象に出会うと「どちらかの解析が間違っているのでは」と疑いたくなりますが、どちらも正しく、単に推定している対象（エスティマンド）が違うだけです。
・GEEは周辺（marginal）効果＝集団平均の効果を推定しています。臨床的な問いに翻訳すると「この治療を集団全体に適用したら、集団の中で状態良好となる割合は平均的にどれだけ変わるか」。公衆衛生的・政策的な問い、そして多くの臨床試験の主要な問いはこちらです。
・GLMMは条件付き（conditional）効果＝個人内の効果を推定しています。翻訳すると「同じ患者が治療を受けたら、その人のオッズはどれだけ変わるか」。ランダム切片で表現された「その患者固有のかかりやすさ」を固定したうえでの効果です。
そして重要なのは、ロジットのような非線形リンクを使うかぎり、周辺効果は条件付き効果より必ず0（＝OR 1）の方向に縮む（attenuation、希釈）という性質があることです。個人ごとに異なる切片を持つロジスティック曲線を集団全体で平均すると、平均された曲線は個々の曲線より傾きがゆるやかになるためです。しかもこの縮みは、ランダム効果のばらつきが大きいほど大きくなります。今回はランダム切片の標準偏差が 1.99（分散 3.96）とかなり大きく、患者間の「状態良好になりやすさ」の差が非常に大きいデータでした。1.26 と 2.12 という大きな乖離は、まさにこの大きな個人差が生んだものです。
逆に言えば、線形モデル（恒等リンク）では周辺効果と条件付き効果が一致します。平均をとる操作と線形変換は順番を入れ替えられるからです。つまり、連続量アウトカムのMMRMを使っているかぎり、この区別は表面化しません。生物統計の実務でこの論点が突然重くなるのは、二値・カウントといった非線形リンクのアウトカムを反復測定で扱うときだけだ、と覚えておくと整理しやすいはずです。
以上を踏まえ、3つの手法の位置づけを一覧に整理します。
| 手法 | 推定する効果 | 主な用途 | 欠測への仮定 | Rの実装 |
|---|---|---|---|---|
| GEE | 周辺効果（集団平均） | 二値・カウントの反復測定データ。集団としての平均的な治療効果を示したいとき | MCAR（完全にランダムな欠測）を要求 | geepack::geeglm |
| GLMM | 条件付き効果（個人内） | 個人差そのものに関心があるとき。変量効果の分散を推定したいとき。施設差・患者差を明示的に分解したいとき | MAR（ランダムな欠測）で妥当 | lme4::glmer |
| MMRM | 連続量のため周辺効果と条件付き効果が一致 | 連続量の反復測定データ。規制当局にもっとも広く受け入れられている連続量の主解析 | MAR（ランダムな欠測）で妥当 | mmrm::mmrm |
実務的な使い分けの結論はシンプルです。連続量の反復測定であればMMRMが第一選択で、規制当局への説明コストがもっとも低い選択肢になります（実装はmmrmをRで実装するで詳しく扱っています）。二値・カウントの反復測定で「集団としての治療効果」を報告したいならGEE、個人差の大きさそのものが関心事ならGLMM、という整理になります。
表の「欠測への仮定」の列は、GEEを臨床試験で使うときの最大の弱点です。通常のGEEはMCAR（Missing Completely At Random、完全にランダムな欠測）を仮定しており、脱落がアウトカムそのものに関連している状況、すなわちMAR（Missing At Random）の状況では推定が偏ります。
これは机上の話ではありません。臨床試験では「効果が乏しい患者ほど途中で脱落する」「有害事象が出た患者ほど来院しなくなる」といった、アウトカムと関連した脱落が起きるのがむしろ普通です。この状況で通常のGEEをそのまま主要解析に据えると、治療効果を過大にも過小にも見せてしまう恐れがあります。
対処としては、脱落確率の逆数で観測を重み付けする重み付きGEE（IPW-GEE, inverse probability weighted GEE）や、多重代入法（multiple imputation）で欠測を補完してからGEEを当てはめる方法が使われます。欠測の扱いという論点全般については、MMRM（反復測定混合モデル）と多重代入法の組み合わせ解析で考え方を整理していますので、そちらもあわせて参照してください。
実務でのポイント
・クラスタ数が少ないときは頑健分散を過信しない：頑健分散は「クラスタ数が十分に多い」という漸近論に支えられています。目安として40クラスタ未満では標準誤差が過小推定になりやすく、その結果として第一種の過誤が名目の5%を超えてしまいます。少数施設のクラスタランダム化試験などではKauermann–Carroll型やMancl–DeRouen型といったスモールサンプル補正の適用を検討してください。
・解析計画書（SAP）で作業相関構造を事前規定する：作業相関構造は「複数試してあてはまりの良かったものを選ぶ」ものではありません。データを見てから選ぶと多重性の問題が生じます。試験の来院間隔やアウトカムの性質から事前に決め、SAPに明記し、他の構造は感度分析として位置づけるのが実務の作法です。
・「集団平均の効果」を報告していることを明記する：GEEの結果を報告するときは、それが周辺（集団平均）効果であることを本文や脚注で明示します。同じデータのGLMMと数値が違っても矛盾ではありませんが、その説明を用意しておかないと査読やレビューで必ず指摘されます。
・欠測がMARの疑いがあるなら通常のGEEをそのまま主解析にしない：脱落パターンをまず記述統計で確認し、アウトカムと関連していそうなら重み付きGEEや多重代入法との併用を検討します。「とりあえずGEE」で押し切らないことが重要です。
・クラスタIDの作り方を必ず検算する：本記事で interaction(center, id) を作ったように、多施設試験の患者IDは施設内で振り直されていることが少なくありません。出力の Number of clusters が想定している被験者数と一致しているかを、解析のたびに必ず目視で確認してください。ここを間違えると、モデル式が正しくても結果全体が意味を失います。
📚 この記事をより深く理解するための参考書籍
関連記事
GEEはGLMの拡張として位置づけられる手法です。土台となる考え方や、二値以外のアウトカムへの応用については以下の記事もあわせてご覧ください。
- 一般化線形モデル(GLM)の“裏側” ― 指数型分布族とリンク関数から、GEEの土台となるGLMの構造を解説しています
- ポアソン回帰とは ― 件数・発生率データの解析。GEEをカウントアウトカムに使う前提となる知識です
- 負の二項回帰と再発イベント率解析 ― 過分散のあるカウントデータの扱い。喘息・COPD増悪のような再発イベントの解析に直結します
まとめ
本記事では、GEE（一般化推定方程式）を「周辺モデル」という考え方から出発して整理してきました。個体差を変量効果として明示的にモデル化するのではなく、相関構造は作業仮説として扱い、平均構造だけを正しく推定するというGEEの設計思想は、一度腹に落ちるとその後の判断が一気に楽になります。作業相関構造の指定が多少ずれていても、頑健分散（サンドイッチ分散）によって標準誤差の妥当性が守られる。この「一致性の担保の仕方」こそがGEEの本質でした。
実装面では、geepack::geeglm を使えば id と corstr を指定するだけでGEEが動くこと、そして通常の glm と比べたときに、被験者間で一定の共変量では標準誤差が広がり、被験者内で変化する共変量ではむしろ狭くなるという非対称な変化が起きることを、実データの出力で確認しました。QICによる構造の比較も含め、Rでの一通りの手順は追えたはずです。そして本記事最大の学びどころは、同じデータ・同じ共変量でありながら、GEEのオッズ比 3.53 とGLMMのオッズ比 8.37 が2倍以上も食い違ったところにあります。これは誤りではなく、周辺効果と条件付き効果という異なるエスティマンドを推定していることの必然的な帰結でした。
臨床試験の現場では、「この治療を集団に適用したら平均的にどうなるか」という問いに答えることが求められる場面が数多くあります。二値やカウントのアウトカムを反復測定で扱うとき、その問いに正面から答えられるのがGEEです。同時に、MCAR仮定という弱点を理解し、脱落がアウトカムに関連しそうな場面では重み付きGEEや多重代入法を検討する。手法の長所と限界の両方を語れることは、生物統計家として議論をリードするうえで大きな強みになります。
