一般化推定方程式(GEE)の基礎と臨床研究での応用

臨床研究や疫学研究において、繰り返し測定されるデータやクラスター化されたデータは頻繁に登場する。このようなデータでは、同じ被験者からの測定値や同じクラスター内の観測値には相関があるため、従来の独立性を仮定する統計手法では適切な解析ができない。一般化推定方程式（Generalized Estimating Equations: GEE）は、このような相関のあるデータに対して、ロバストな推定を可能にする強力な統計的手法である。本稿では、GEEの概要から、類似する手法である一般化線形混合モデル（GLMM）との違い、具体的な適用場面、4つの異なるデータタイプに応じた臨床研究での応用例、そしてそれぞれのRを用いた計算例までを網羅的に解説する。

一般化推定方程式の概要

一般化推定方程式（GEE）は、対象内で相関のあるデータ（例えば、同じ患者から繰り返し測定されたデータや、家族内で複数のメンバーから得られたデータなど）の解析に用いられる統計的手法である。GEEは、平均構造（目的変数の平均が予測変数によってどのように変化するか）をモデル化することに焦点を当て、相関構造を「ニューサンス・パラメータ」（後述）として扱うことで、相関構造が誤って指定された場合でも回帰係数の推定値がロバスト（頑健）であることが特徴である。

GEEは、一般化線形モデル（GLM）の枠組みを拡張したものであり、応答変数の分布が正規分布に限らず、二項分布、ポアソン分布、順序ロジットなど、指数型分布族に属する様々な分布に対応できる。推定は、各対象内で（本当は独立ではないが、あえて）独立な観測値が存在するという仮定の下で得られる標準誤差を、サンドイッチ推定（ロバスト標準誤差）によって調整することで行われる。これにより、相関の具体的な構造を正確に特定する必要がなく、解析の柔軟性が高まる。

ニューサンス・パラメータとしての相関構造

GEEにおける「ニューサンス・パラメータ（Nuisance Parameter）」（局外パラメータとも言われる）としての相関構造の扱いは、この手法の重要な特徴であり、強みの一つである。

統計モデリングにおいて、関心のある主要なパラメータ（例えば、治療効果を示す回帰係数）とは別に、モデルの適合に必要ではあるものの、直接的な研究目的の対象ではないパラメータを「ニューサンス・パラメータ」と呼ぶ。GEEでは、データ内の相関構造（例：繰り返し測定間の相関のパターン）がこれに該当する。

GEEの設計思想は、「平均構造を正しく指定すれば、相関構造が完璧に指定されていなくても、回帰係数の推定値は一致性を保つ（正しい値に近づく）」という点にある。これは、GEEが、各観測単位（例：患者）内の相関を明示的にモデル化するのではなく、その存在を「ロバスト標準誤差」という形で補正することで考慮する。

具体的には、GEEは以下のステップで推定を行う。

1. まず、各観測単位内のデータが独立であるという「作業独立性」を仮定して、平均構造の回帰係数を暫定的に推定する。
2. この暫定的な回帰係数を用いて、観測されたデータ間の相関を推定する（例えば、「交換可能」や「AR(1)」といった仮定に基づき、クラスター内の相関を計算する）。
3. この推定された相関構造を用いて、回帰係数を再推定する。このプロセスを収束するまで繰り返す。
4. 最終的な回帰係数の推定値が得られたら、その標準誤差を「サンドイッチ推定量（またはロバスト標準誤差）」という特殊な方法で算出する。このサンドイッチ推定量によって、相関構造の仮定が多少間違っていても、標準誤差の推定を頑健にすることができる。

このアプローチにより、分析者は相関構造の複雑な詳細に過度にこだわる必要がなくなる。例えば、繰り返し測定データにおいて、測定時点が不規則であったり、欠損データが多かったりする場合でも、GEEは比較的容易に適用できる。この柔軟性とロバスト性が、GEEが臨床研究で用いられる理由となっている。

一般化線形混合モデル (GLMM) との違い

GEEと一般化線形混合モデル（GLMM）は、いずれも相関のあるデータに対応できる強力な手法であるが、そのアプローチと目的には重要な違いがある。

特徴	一般化推定方程式 (GEE)	一般化線形混合モデル (GLMM)
アプローチ	母集団平均に対する回帰係数の推定（集団平均効果）	個体差を考慮した回帰係数の推定（個体レベルでの効果）
推定の焦点	平均構造の推定に主眼。相関構造はニューサンス・パラメータとして扱う。	平均構造と分散共分散構造の両方を明示的にモデル化する。
相関構造	誤って指定しても推定値はロバスト。	明示的に指定する必要があり、誤ると推定に偏りが生じる可能性がある。
欠損データ	欠損メカニズムがMCAR (Missing Completely at Random) または MAR (Missing at Random) の場合に対応可能。MCARが理想的。	MAR (Missing at Random) の欠損データに対応しやすい。
仮定	観測値の独立性に対する仮定は不要。ロバスト標準誤差を用いる。	ランダム効果が正規分布に従うという仮定などが必要。
主な用途	集団レベルの効果を推定したい場合。	個体レベルの変動や個体間の異質性を考慮したい場合。

簡単に言えば、GEEは「集団全体としてどうか」という問いに答えるのに適しており、GLMMは「個々の個体がどのように異なるか」という問いに答えるのに適していると言える。

一般化推定方程式の使い所

GEEは以下のような状況で特に有用である。

• 繰り返し測定データ: 同じ被験者から複数回にわたって測定されたデータ（例：薬物投与後の血圧の経時変化、治療効果の追跡）。
• クラスター化データ: 家族、学校、病院など、特定のグループ内に属する個人からのデータ（例：同じ家族内のメンバーの健康状態、同じ学校の生徒の学習成績）。
• 縦断研究: 長期間にわたって追跡調査される研究で、被験者内の相関を考慮する必要がある場合。
• 集団レベルの効果推定: 個々の変動よりも、介入や因子の集団全体への平均的な効果に関心がある場合。
• 非正規分布の応答変数: 二項データ（成功/失敗）、カウントデータ（イベント発生数）、順序尺度データなど、応答変数が正規分布に従わない場合。

臨床研究の具体例とRによる計算例

ここでは、4つの異なるデータタイプ（連続データ、2項アウトカム、順序カテゴリ、カウントデータ）について、GEEの具体的な適用例とRによる計算例を示す。

連続データ：高血圧治療薬の効果検証

• 目的: 新しい降圧剤の血圧降下効果を、従来薬と比較して評価する。
• データ: 複数の患者を新薬群と従来薬群に割り付け、治療開始前、1ヶ月後、3ヶ月後、6ヶ月時点での収縮期血圧を測定。
• GEEの適用: 患者内の血圧測定値の相関（繰り返し測定）を考慮し、新薬の平均的な効果を従来薬と比較して推定する。応答変数は連続変数（正規分布を仮定）。

R による計算例：

# パッケージのインストールと読み込み
# install.packages("geepack")
library(geepack)

# 架空のデータ作成 (連続データ)
set.seed(123)
id <- rep(1:100, each = 4) # 患者ID (100人)
time <- rep(c(0, 1, 3, 6), times = 100) # 測定時点 (0ヶ月, 1ヶ月, 3ヶ月, 6ヶ月)
# 治療グループ (0:従来群, 1:新薬群) - ランダムに割り付け
trt <- rep(sample(c(0, 1), 100, replace = TRUE), each = 4)

# 収縮期血圧の生成 (ベースラインは高め、時間経過と治療効果で低下)
# 従来群: ベースライン140 mmHg, 時間経過で緩やかに低下
# 新薬群: ベースライン140 mmHg, 時間経過でより大きく低下
sbp <- numeric(length(id))
for (i in 1:length(id)) {
    base_sbp <- rnorm(1, 140, 5) # 個体ごとのベースラインのばらつき
    if (trt[i] == 0) { # 従来群
        sbp[i] <- base_sbp - 2 * time[i] + rnorm(1, 0, 3)
    } else { # 新薬群
        sbp[i] <- base_sbp - 5 * time[i] + rnorm(1, 0, 3)
    }
}
# 血圧が極端に低くならないように調整
sbp[sbp < 100] <- 100 + rnorm(sum(sbp < 100), 0, 2)

data_sbp <- data.frame(id, time, trt, round(sbp))
head(data_sbp)

# GEEモデルの実行 (連続データ)
# family = gaussian: 応答変数が連続であるため
# link = "identity": 恒等リンク関数
model_sbp_gee <- geeglm(sbp ~ trt + time + trt:time,
    id = id,
    data = data_sbp,
    family = gaussian(link = "identity"),
    corstr = "exchangeable"
) # 交換可能な相関構造を仮定

# 結果の表示
summary(model_sbp_gee)

実行結果：

> # 結果の表示
> summary(model_sbp_gee)

Call:
geeglm(formula = sbp ~ trt + time + trt:time, family = gaussian(link = "identity"), 
    data = data_sbp, id = id, corstr = "exchangeable")

 Coefficients:
            Estimate  Std.err     Wald Pr(>|W|)    
(Intercept) 140.5775   0.5688 61089.83   <2e-16 ***
trt          -0.8562   0.8976     0.91     0.34    
time         -2.1036   0.1609   170.90   <2e-16 ***
trt:time     -2.7607   0.2770    99.36   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation structure = exchangeable 
Estimated Scale Parameters:

            Estimate Std.err
(Intercept)    30.98   2.114
  Link = identity 

Estimated Correlation Parameters:
      Estimate Std.err
alpha -0.02664 0.04149
Number of clusters:   100  Maximum cluster size: 4 
> 

結果の見方としては、trt（治療グループ）、time（測定時点）、trt:time（治療グループと時間の交互作用）の係数を確認する。特にtrt:timeの係数が統計的に有意であれば、新薬群と従来群で血圧の経時変化に有意な差があることを示唆している。今回の結果では、交互作用が有意であるため、新薬群と従来群で血圧の経時変化に差異があることが示唆されているが、主効果自体は、有意ではなく、新薬群の効果が示唆されたわけではない。

2項アウトカム：新規治療法による疾患寛解率の比較

• 目的: 新規治療法が、ある疾患の寛解（病状が一時的または永続的に消失すること）達成率を、従来治療法と比較して向上させるかを評価する。
• データ: 患者を新規治療群と従来治療群に割り付け、治療開始から1ヶ月、3ヶ月、6ヶ月後に寛解の有無（0:なし, 1:あり）を評価。
• GEEの適用: 各患者の寛解状況の経時変化における相関を考慮し、新規治療法と従来治療法の寛解達成率の群間差を評価する。応答変数は二値。

R による計算例：

# パッケージのインストールと読み込み
# install.packages("geepack")
library(geepack)

# 架空のデータ作成 (2項アウトカム)
set.seed(456)
id_binary <- rep(1:80, each = 3) # 患者ID (80人)
time_binary <- rep(c(1, 3, 6), times = 80) # 測定時点 (1ヶ月, 3ヶ月, 6ヶ月)
# 治療グループ (0:従来群, 1:新規治療群)
trt_binary <- rep(sample(c(0, 1), 80, replace = TRUE), each = 3)

# 寛解の有無 (0:なし, 1:あり) の生成
# 新規治療群は寛解率が高く、時間経過で寛解しやすくなる
remission <- numeric(length(id_binary))
for (i in 1:length(id_binary)) {
    prob_remission <- 0.2 # ベースラインの確率
    if (trt_binary[i] == 1) { # 新規治療群
        prob_remission <- prob_remission + 0.1 * time_binary[i] + 0.2 # 時間と治療効果
    } else { # 従来群
        prob_remission <- prob_remission + 0.05 * time_binary[i]
    }
    remission[i] <- rbinom(1, 1, min(prob_remission, 0.9)) # 確率が1を超えないように
}

data_remission <- data.frame(id = id_binary, time = time_binary, trt = trt_binary, remission)
head(data_remission)

# GEEモデルの実行 (2項アウトカム)
# family = binomial: 応答変数が二値であるため
# link = "logit": ロジットリンク関数
model_remission_gee <- geeglm(remission ~ trt * as.factor(time),
    id = id,
    data = data_remission,
    family = binomial(link = "logit"),
    corstr = "exchangeable"
)

# 結果の表示
summary(model_remission_gee)


OR <- exp(coef(model_remission_gee))
CI_lower <- exp(coef(model_remission_gee) - 1.96 * sqrt(diag(vcov(model_remission_gee)))) # vcov ロバスト分散共分散行列
CI_upper <- exp(coef(model_remission_gee) + 1.96 * sqrt(diag(vcov(model_remission_gee)))) # vcov ロバスト分散共分散行列
p_value <- ifelse(summary(model_remission_gee)$coefficients[, 4] < 0.001, "<0.001", round(summary(model_remission_gee)$coefficients[, 4], 3))

odds_ratios <- data.frame(
    OR = round(OR, 2),
    CI_lower = round(CI_lower, 2),
    CI_upper = round(CI_upper, 2),
    p_value = p_value
)

print("オッズ比とその95%信頼区間:")
print(odds_ratios)

実行結果：

[1] "オッズ比とその95%信頼区間:"
> print(odds_ratios)
                       OR CI_lower CI_upper p_value
(Intercept)          0.33     0.16     0.68   0.003
trt                  3.67     1.42     9.47   0.007
as.factor(time)3     2.00     0.77     5.20   0.155
as.factor(time)6     3.32     1.34     8.20   0.009
trt:as.factor(time)3 1.23     0.34     4.45   0.755
trt:as.factor(time)6 2.22     0.43    11.36   0.338

結果の解釈としては、trtの係数が統計的に有意であれば、新規治療群と従来群の間で寛解率に有意な差があることを示唆している。また、交互作用項 trt:as.factor が有意であれば、時間的な変化が異なることを示している。今回は、交互作用は有意ではなく、時間によって群間差が異なるとは言えない。trt は オッズ比 3.67 と計算され、統計学的有意に異なる（P = 0.007）という結果であった。

順序カテゴリ：疼痛スコアの経時的変化の比較

• 目的: 新しい治療法が、呼吸器疾患患者の全般的症状レベルを、プラセボと比較して効果的に改善するかを評価する。
• データ: 患者を新治療群とプラセボ群に割り付け、全4回の来院時（y1, y2, y3, y4）に全般的症状レベルを3段階（1 = poor（悪い）、2 = good（良い）、3 = excellent（非常に良い）で評価。
• GEEの適用: 各患者の全般性症状レベルの経時変化における相関を考慮し、両群間の症状改善効果を比較する。応答変数は順序カテゴリ。

Rによる計算例:

library(geepack)
data(respdis)
head(respdis)

# データの整形(ワイド形式からロング形式に変換)
resp.l <- reshape(respdis,
    varying = list(c("y1", "y2", "y3", "y4")),
    v.names = "resp", direction = "long"
)
head(resp.l)

# データの並べ替え(idとtimeでソート)
resp.l <- resp.l[order(resp.l$id, resp.l$time), ]
head(resp.l)

# データの確認(time, resp, trtのクロス表)
table(resp.l$time, resp.l$resp, resp.l$trt)

# モデルの推定(時間と処置の交互作用を考慮)
fit <- ordgee(ordered(resp) ~ trt * time, id = id, data = resp.l, corstr = "exchangeable")
summary(fit)

# オッズ比の計算
OR <- exp(fit$beta)
CI_lower <- exp(fit$beta - 1.96 * sqrt(diag(fit$vbeta))) # vbeta ロバスト分散共分散行列
CI_upper <- exp(fit$beta + 1.96 * sqrt(diag(fit$vbeta))) # vbeta ロバスト分散共分散行列
P <- 1 - pchisq((fit$beta / sqrt(diag(fit$vbeta)))^2, df = 1)
p_value <- ifelse(P < 0.001, "<0.001", round(P, 3))

odds_ratios <- data.frame(
    OR = round(OR, 2),
    CI_lower = round(CI_lower, 2),
    CI_upper = round(CI_upper, 2),
    p_value = p_value
)

print("オッズ比とその95%信頼区間:")
print(odds_ratios)

実行結果：

[1] "オッズ比とその95%信頼区間:"
> print(odds_ratios)
           OR CI_lower CI_upper p_value
Inter:1  3.61     2.10     6.21  <0.001
Inter:2  0.37     0.21     0.65  <0.001
trt      3.59     1.42     9.07   0.007
time     0.95     0.81     1.11   0.508
trt:time 0.92     0.70     1.21   0.541

まず最初に目に入る Inter:1, Inter:2 であるが、これは、アウトカムの 3つのカテゴリの1番目のカットポイント（1を超える、つまり2または3になるポイント）と2番めのカットポイント（2を超える、つまりカテゴリ3になるポイント）を示している。これは表示されるのだが、このカットポイントは提示が必要な値ではなく、解釈も不要である。1，2，3のカットポイントの幅は「1」ではないことが特徴で、その点がノンパラメトリック的（つまり順序に過ぎない）なんだと理解すれば良い。

trt, time, trt:time が見るべき行である。trt が1より大きく統計学的に有意（P < 0.05）であるため、新治療がプラセボに比べて、症状を改善することを示唆している。time, trt:time は OR が 1 に近く、統計学的有意でもないため、経時的な効果や治療法と時間の交互作用は明らかではないと解釈できる。

カウントデータ：感染症の発生件数比較

• 目的: 新しい予防介入プログラムが、ある集団における感染症の年間発生件数を、従来の介入と比較して削減するかを評価する。
• データ: 複数の対象集団（例：地域コミュニティ、施設など）を新予防プログラム導入群と従来予防策実施群に割り付け。各集団について、ベースライン（介入前1年間）、追跡開始1年後、追跡開始2年後の年間感染症発生件数を記録。
• GEEの適用: 各集団内の経時的な感染症発生件数の相関を考慮し、予防プログラムの全体的な効果を推定する。応答変数はカウントデータ。コミュニティの大きさもオフセット（後述）で考慮。

R による計算例：

# パッケージの読み込み
library(geepack)

# 架空のデータ作成 (カウントデータ) - population を追加
set.seed(2024)
community_id <- rep(1:50, each = 3) # 集団ID (50コミュニティ)
# 測定時点 (0:ベースライン1年間, 1:追跡1年間, 2:追跡2年間)
measurement_period <- rep(c(0, 1, 2), times = 50)
# 予防プログラム (0:従来予防策, 1:新予防プログラム)
program <- rep(sample(c(0, 1), 50, replace = TRUE), each = 3)

# コミュニティの人口 (ベースラインで固定、各コミュニティで異なる)
# 実際のデータでは、各コミュニティの人口を period ごとに持つ場合もありますが、
# この例ではシンプルにコミュニティIDごとに固定とします。
population_base <- rep(sample(5000:20000, 50, replace = TRUE), each = 3)

# 各期間における感染症発生件数の生成
# population の影響を考慮し、感染率を基準に生成
infection_count <- numeric(length(community_id))
for (i in 1:length(community_id)) {
    # 人口あたりのベースライン感染率 (例: 1000人あたり15件)
    # rate_per_1000 = 15
    # base_expected_count = (population_base[i] / 1000) * rate_per_1000
    base_expected_count <- (population_base[i] / 1000) * 15 # ベースラインで人口1000人あたり15件

    if (measurement_period[i] == 0) { # ベースライン
        # ポアソン分布から生成、期待値はベースラインの人口に応じる
        infection_count[i] <- rpois(1, base_expected_count)
    } else if (program[i] == 0) { # 従来予防策群 (追跡期間)
        # 緩やかに減少 (例: 1000人あたり年0.5件減少)
        infection_count[i] <- rpois(1, max(0.1, base_expected_count - (measurement_period[i] * 0.5 * (population_base[i] / 1000))))
    } else { # 新予防プログラム群 (追跡期間)
        # より大きく減少 (例: 1000人あたり年2件減少)
        infection_count[i] <- rpois(1, max(0.1, base_expected_count - (measurement_period[i] * 2 * (population_base[i] / 1000))))
    }
    infection_count[i] <- max(0, infection_count[i]) # 負の数にならないように
}

data_infection_count <- data.frame(
    id = community_id,
    period = measurement_period,
    program = program,
    population = population_base, # 人口を追加
    infection_count
)
head(data_infection_count)

# GEEモデルの実行 (カウントデータ) - オフセットを追加
# family = poisson: 応答変数がカウントデータであるため
# link = "log": 対数リンク関数
# offset = log(population): コミュニティの人口の対数をオフセットとして指定
model_infection_gee_offset <- geeglm(infection_count ~ program + as.factor(period) + program:as.factor(period),
    id = id,
    data = data_infection_count,
    family = poisson(link = "log"),
    corstr = "exchangeable",
    offset = log(population)
) # オフセットの指定

# 結果の表示
summary(model_infection_gee_offset)

RR <- exp(coef(model_infection_gee_offset))
CI_lower <- exp(coef(model_infection_gee_offset) - 1.96 * sqrt(diag(vcov(model_infection_gee_offset)))) # vcov ロバスト分散共分散行列
CI_upper <- exp(coef(model_infection_gee_offset) + 1.96 * sqrt(diag(vcov(model_infection_gee_offset)))) # vcov ロバスト分散共分散行列
p_value <- ifelse(summary(model_infection_gee_offset)$coefficients[, 4] < 0.001, "<0.001", round(summary(model_infection_gee_offset)$coefficients[, 4], 3))

rate_ratios <- data.frame(
    RR = round(RR, 2),
    CI_lower = round(CI_lower, 2),
    CI_upper = round(CI_upper, 2),
    p_value = p_value
)

print("レート比とその95%信頼区間:")
print(rate_ratios)

実行結果：

> # 結果の表示
> summary(model_infection_gee_offset)

Call:
geeglm(formula = infection_count ~ program + as.factor(period) + 
    program:as.factor(period), family = poisson(link = "log"), 
    data = data_infection_count, offset = log(population), id = id, 
    corstr = "exchangeable")

 Coefficients:
                           Estimate  Std.err     Wald Pr(>|W|)    
(Intercept)                -4.18028  0.01472 80688.53  < 2e-16 ***
program                     0.00188  0.01945     0.01  0.92293    
as.factor(period)1         -0.04594  0.02427     3.58  0.05831 .  
as.factor(period)2         -0.08454  0.01449    34.05  5.4e-09 ***
program:as.factor(period)1 -0.12009  0.03280    13.41  0.00025 ***
program:as.factor(period)2 -0.28120  0.02262   154.52  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation structure = exchangeable 
Estimated Scale Parameters:

            Estimate Std.err
(Intercept)     1.04   0.106
  Link = identity 

Estimated Correlation Parameters:
      Estimate Std.err
alpha    0.101    0.09
Number of clusters:   50  Maximum cluster size: 3 
> 

結果の解釈:

programの係数は、ベースライン時点での両群間の平均的な発生件数の差（対数レート比）を示している。as.factor(period)1やas.factor(period)2の係数は、従来群におけるベースラインからの各期間への変化を示している。最も重要なのは交互作用項 program:as.factor(period)1 と program:as.factor(period)2 である。これらの係数が統計的に有意であれば、新予防プログラムが従来予防策と比較して、各追跡期間において感染症発生件数の変化に有意な差をもたらしたことを示唆している。

populationは、各community_idに対応する架空のpopulation（コミュニティの大きさ）を表していて、geeglm関数の引数にoffset = log(population)を追加して、オフセットを設定する。ポアソン回帰では対数リンク関数を使うため、オフセットも対数変換したものを指定する。オフセット項としてlog(population)を含めることで、モデルの係数は人口あたりの感染症発生率（レート）に対する効果として解釈される。

[1] "レート比とその95%信頼区間:"
> print(rate_ratios)
                             RR CI_lower CI_upper p_value
(Intercept)                0.02     0.01     0.02  <0.001
program                    1.00     0.96     1.04   0.923
as.factor(period)1         0.96     0.91     1.00   0.058
as.factor(period)2         0.92     0.89     0.95  <0.001
program:as.factor(period)1 0.89     0.83     0.95  <0.001
program:as.factor(period)2 0.75     0.72     0.79  <0.001

係数は対数レート比で計算されているので、exp()関数でレート比に変換すると直感的な解釈が可能になる。レート比となると、新予防プログラム群の発生率が従来群の発生率の何倍になるかという解釈が可能である。今回は、program のレート比を見ると 1 であるが、交互作用の係数は1より下回っており、統計学的有意なので、従来群と比較して新予防群では、交互作用ありで減少していることが伺える。

オフセットとは？

オフセット（offset）は、GLMMのような一般化線形モデルや GEE において、モデルの線形予測子に既知の固定された値を追加する項のことである。この値は、通常、目的変数の期待値に影響を与えるエクスポージャー（露出）や背景の大きさを表すものである。

上記の例で「コミュニティの大きさ（例：各コミュニティの人口）」が異なる場合、感染症発生件数はコミュニティの規模に依存すると考えられる。この場合、コミュニティの大きさをオフセットとして加えることで、人口1000人あたりの感染症発生率などのレートを比較できるようになる。

上記架空のデータには、populationというコミュニティ規模を示す変数があり、オフセットとして使用した例を示した。

臨床研究におけるGLMM優先の理由（推測）

GEEが非常に有用であるにもかかわらず、臨床研究でGLMM（一般化線形混合モデル）の方が多く見られる。いくつか理由が考えられるが、主に以下の点が背景にあると推測される。

個体レベルの異質性（Subject-Specific Effects）への関心

臨床研究、特に薬剤開発や個別化医療の文脈では、個々の患者の反応のばらつきや、治療効果が患者によってどのように異なるかに関心が高いことが多い。GLMMは、ランダム効果を導入することで、このような個体間の異質性を直接モデル化し、患者ごとの予測や効果を推定することができる。

一方、GEEは集団平均効果（Population-Average Effects）の推定に焦点を当てている。これは、ある治療が「平均的に」どの程度効果があるかを知るには適しているが、個々の患者の治療計画を立てる上では、個別の反応を考慮できるGLMMの方がより情報量が多いと感じられるのかもしれない。

なぜGEEよりGLMMか: 多くの臨床試験は、治療が「誰に」どの程度効くかという問いに関心がある。GLMMは、ランダム効果を通じて、この「誰に」という要素を説明できる。例えば、ある患者はベースラインの特性が似ていても、別の患者より薬が効きにくい、あるいは効きやすいといった個体差をモデルに組み込めるのはGLMMの強みである。

欠損データへの対応能力

臨床研究では、患者のドロップアウトや測定忘れなどにより、欠損データが頻繁に発生する。

• GLMM: MAR (Missing At Random) の欠損メカニズムに対応する能力が高いとされている。これは、欠損が観測されたデータ（例：これまでの測定値やベースラインの特性）に依存するが、欠損した値そのものには依存しない場合を指す。GLMMは、欠損データを尤度ベースで自然に扱えるため、追加の調整なしにMARの欠損に対応しやすい。
• GEE: GEEは、基本的には欠損がMCAR (Missing Completely At Random) である場合に最適に機能する。MCARは、欠損が完全にランダムで、観測されたデータも欠損したデータも、いずれにも依存しない状況を指す。もし欠損がMARの場合、GEEは単純に欠損した観測を削除（完全ケース解析）するとバイアスが生じる可能性がある。GEEでMARの欠損に対応するには、逆確率重み付け（Inverse Probability Weighting: IPW）などのより複雑な手法と組み合わせる必要があり、解析のハードルが上がる。

なぜGEEよりGLMMか: 臨床研究における欠損データはMARであることが多く、GLMMの方がデフォルトでよりロバストに対応できるため、解析の手間やバイアスの懸念が少ないと認識されているからかもしれない。

モデルの解釈と直感的な理解

GLMMは、固定効果と変量効果という明確な枠組みでデータ構造を表現するため、統計専門家でなくとも、その概念が比較的理解しやすいと感じられるのかもしれない。変量効果は、「個体差がどれくらいあるか」を分散成分として直接的に推定するため、データの背後にあるメカニズムをより詳細に記述していると理解できる。

GEEの「ニューサンス・パラメータとしての相関構造」という考え方は、数理統計学的な背景がないと直感的に理解しにくい可能性がある。相関構造を間違えても係数がロバストである、という利点は強力であるが、その「仕組み」が分かりにくいと感じられるのかもしれない。

なぜGEEよりGLMMか: 「ランダム効果」という概念は、患者間のばらつきを表現するのに適しており、臨床家や研究者にとって直感的に受け入れられやすい。一方、GEEの頑健性の根拠である「サンドイッチ推定量」の概念は、より専門的な知識を要する。

統計ソフトウェアの普及と使いやすさ

GLMMは、SASのPROC MIXEDやPROC GLIMMIX、Rのlme4パッケージなど、主要な統計ソフトウェアで非常に使いやすく実装されており、多くのチュートリアルや解説資料が存在する。これにより、研究者がGLMMを学習し、適用する敷居が低いという側面があるだろう。

GEEもRのgeepackやgeeパッケージなどで利用できるが、解説資料は比較的少ない。

なぜGEEよりGLMMか: ソフトウェアの普及度や、一般的な統計解析のトレーニングでGLMMがより多く扱われる傾向にあるため、研究者がGLMMを選択しやすくなっている可能性がある。

GEEかGLMMか

GEEのロバスト性は非常に魅力的であるが、臨床研究の多くが「個体レベルの治療効果や反応のばらつき」に関心を持ち、欠損データへの対応やモデルの直感的な解釈、そして既存の統計ソフトウェアの使いやすさといった実用的な側面から、GLMMが選ばれることが多いと考えられる。

しかし、これはGEEが劣っているという意味ではない。集団全体の効果を頑健に評価したい場合や、相関構造を厳密に特定することが難しい場合には、GEEは依然として優れた選択肢である。研究の目的とデータの特性に応じて、適切な統計手法を選択することが最も重要である。

まとめ

一般化推定方程式（GEE）は、繰り返し測定データやクラスター化データなど、相関のあるデータに対してロバストな回帰分析を可能にする重要な統計的手法である。GLMMとは異なり、個体レベルの変動ではなく集団平均効果の推定に焦点を当て、相関構造の誤指定に頑健であるという利点がある。この頑健性は、相関構造を「ニューサンス・パラメータ」として扱い、ロバスト標準誤差を用いることで実現される。

本稿では、GEEの概要からGLMMとの違い、そして連続データ、2項アウトカム、順序カテゴリ、カウントデータといった異なるタイプの応答変数に対する臨床研究の具体例とRによる計算例を示した。これらの例から、GEEが様々な種類の相関データを適切に解析するための強力なツールであることが理解いただければ幸いである。研究デザインと目的に応じてGEEの適切な適用を検討することが、より信頼性の高い研究結果を得る上で不可欠である。